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Chapter 1 


Introduction 


1.1 Problem Statement 

The objective of this study is to numerically examine the interaction of a premixed laminar 
flame with its self-induced boundary layer in a channel. As a flame propagates down a channel, 
it creates pressure waves ahead of it . 1 If the flame is accelerating, the flame generates compression 
waves and if the flame is decelerating, expansion waves are formed. These waves cause motion of 
the gas ahead of the flame. If the channel has an open end toward which the flame propagates, the 
gas set in motion by the flame can exit the channel. This self-induced fluid motion ahead of the 
flame causes a boundary layer to form along the channel wall. The flame then propagates into this 
boundary layer created by its own motion. 

A flame, or a deflagration wave, is a self-sustaining combustion zone moving at subsonic speeds. 
If the combustion zone moves at a supersonic speed, the wave is a detonation. Many factors influence 
the speed of the flame . 2 The effect of the fuel type is important, with different fuels propagating 
at different speeds. Another factor influencing flame speed is the mixture composition, whether the 
mixture is rich, lean, or stoichiometric. A stoichiometric reaction is one where all of the fuel and 
oxidizer is consumed. A rich mixture is one where there is more fuel present than the stoichiometric 
condition, and a lean mixture is one where there is less fuel than the stoichiometric conditions. The 
initial temperature and pressure also affect the flame speed. If the initial mixture is heated, the 
resulting flame speed will be faster. Pressure effects vary with the order of the reaction. The reaction 
order depends on the exponential powers of the reactants in the rate equation. For zero and first 
order reactions, the flame speed increases as the pressure decreases. For second order reactions, the 
flame speed is independent of the pressure, and for reaction orders greater then two, the flame speed 
decreases with decreasing pressure. 
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Two types of flames exist, premixed and diffusion. A premixed flame has the fuel and oxidizer 
mixed at the molecular level before the material is ignited. An example of this type of flame would 
be a Bunsen burner. The fuel and air enter at the bottom of the Bunsen tube and mix before they 
reach the top of the tube where the flame is located. In a diffusion flame, also called a non-premixed 
flame, the fuel and oxidizer are initially separated. The reaction occurs in the zone where the fuel 
and oxidizer come together and mix. A candle is an example of a diffusion flame. The wax of the 
candle is melted, evaporates, diffuses into the air, and bums in the diffusion zone. 

The main purpose of this investigation is to examine the basic flow physics of the formation of 
the boundary layer ahead of the flame and the interaction with the flame as it propagates into the 
boundary layer. Flames in channels with an open end have applications to the flame spread in ducts 
and in mines. It also deals with the flame spread in an aircraft fuel tank. If the tank becomes 
ruptured and ignites, the gases can escape out from the tank through the rupture. The rate of 
the flame propagation through the tank will depend on where the rupture occurs, either ahead of 
the flame or behind it. In the study of the transition from a deflagration wave to a detonation 
wave, localized ignitions have occurred ahead of the flame in the boundary layer . 2-5 The study of 
micro-propulsion, the use of small thrusters on light-weight spacecraft, depends on the combustion 
in tiny chambers and channels to generate the high pressure and high temperature gas for use in the 
nozzle. The effects of the boundary layer on the flame structure is of interest in the micro-channel 
environment of these propulsion devices. 

In this study, the two-dimensional Navier-Stokes equations are used to simulate the propagation 
of a flame in a channel. A single step, two species chemistry model is used which simulates the 
reaction of an acetylene-air mixture. The boundary condition for the temperature at the wall is 
examined for adiabatic and isothermal conditions to study its effect on the flame propagation. The 
flame ignition method is varied from a planar ignition to a spark ignition to examine what effects 
ignition has on the flame structure. 

1.2 Previous Work 

There are three channel configurations that are used for the numerical and experimental work 
dealing with flames in tube. The different configurations depend on the conditions at the ends of the 
tube. Using the notation that the ignition occurs on the left side of the channel, and that the channel 
is described by defining the left and then the right hand sides, the three channel configurations are 
closed-closed, open-closed, and closed-open. 
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1.2.1 Closed - Closed 


Figure LI shows the channel configuration for the closed-closed geometry. The ignition occurs 
at the left closed end and the flame propagates toward the right closed end. A great amount 



Figure 1.1: Closed-Closed Channel Configuration. 


of experimental work 6-12 and numerical work 13 ” 23 has been done on tubes with a closed-closed 
configuration. Figure 1.2 shows the experimental stroboscopic flame record of Ellis 6 for the flame 
propagation in tubes with the same diameter but varying the tube length. For all the tubes, the 
flame was ignited using a spark on the left end. After the sparked ignition, the flame expands 
spherically until it approaches the top and bottom walls. The flame then takes on an elliptical shape 
as the portions that are nearing the top and bottom walls slow down. When the flame reaches these 
walls, the flame is quenched (extinguished) and the flame surface area decreases. This causes the 
flame to slow, since a lower flame area results in less burning, and therefore less hot, expanding gas. 

After the flame slows, one of two possibilities can happen, depending on the length to diameter 
ratio, L/D , of the tube. For tubes with a large L/D ratio, the portions of the flame close to the 
wails begin to move faster than the flame near the centerline. The flame develops a rearward cusp 
near the centerline as the flame at the centerline lags behind the flame closer to the top and bottom 
walls. Salamandra et al. 8 called this a ‘tulip' flame. For tubes with an L/D < 2, no tulip flame 
appears. 6,9 This is seen in figure 1.2 as photograph (d). Photographs (a), (b), and (c) all develop 
the tulip flame. 

The reason for the formation of the tulip is not yet clearly understood. It is known that the 
formation depends on the tube aspect ratio, L/D, the mixture composition, and the initial mixture 
pressure. 9 An earlier explanation by Guenoche 9 theorizes that the tulip is caused by the interaction 
of pressure waves that have reflected off the end wall. During the initial expansion of the flame due 
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Figure 1 .2: Flame propagation in a tube closed at both ends, from Ellis. 6 (a) 19.5 cm, (b) 17 cm, 
(c) 12 cm, (d) 9.5 cm 
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to the spark, the flame emits pressure waves. 1 These waves reflect off of the closed end of the tube 
and intersect the flame. This explanation is dismissed since numerical studies have been performed 
using low Mach number approximations that negate pressure waves and these studies have still 
seen the formation of the tulip flame phenomena. 12, 15,18 Dunn-Rankin et al. 12 determined that 
the presence of vorticity is not required for the tulip formation. It has also been determined that 
viscosity is not required since much numerical work has been done using the inviscid formulation 
and has still generated a tulip flame. 12,16,18,19 The work of Marra and Continillo 23 disagrees with 
this conclusion. They show that when the wall boundary is the slip condition, no tulip forms, but 
if viscosity is included and the condition is the no-slip wall, a tulip forms in the same tube. The 
most recent theory is that the tulip is formed from the Darrieus-Landau 24, 25 instability caused by 
the flow through a curved flame front. 12,18,19,26 

1.2.2 Open - Closed 

Figure 1.3 shows the channel configuration for the open-closed geometry. The ignition is on the 
left hand side and the flame propagates to the right toward the closed end of the tube. The hot 
expanding gases are allowed to escape from the open end of the channel, providing a nearly stagnant 
gas into which the flame moves. 



This open-closed configuration has been studied both experimentally 7, 27-31 and recently numer- 
ically. 16,26,32-34 Mallard and Le Chatelier 27 discovered that for this configuration, the flame travels 
for some distance at a constant speed. Coward and Hartwell 28 found that below a minimum diame- 
ter, the flame speed is not constant. As the diameter is increased, the flame speed becomes constant 
over some distance. As the diameter increases, this constant flame speed also increases. 28 
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Table 1.1: Experimental tubes for the open-closed configuration. 29 


tube number 

i 

2 

3 

4 

Length (cm) 

56 

57 

98 

100 

diameter ^ (cm) 

2.9 

2.9 

4.0 

2.0 

ignition source 

flame 

spark 

spark 

spark 


Figure 1.4, taken from Guenoche and Jouy, 29 shows the flame propagation in various tubes for 
acetylene-air mixtures. Table 1.1 gives the tube parameters for the four tubes used in the experiment. 

In the figure, the flames move from the right to the left. The l h’ and V on the right side of the 
figure refer to the plane that the photographs were taken in, ‘h’ in the horizontal plane and V in 
the vertical. The variables d{ and dj in the original caption refer to the size of the openings on 
the ignition and opposite ends of the tube. The flames are ignited using either a spark or using an 
already developed flame that enters the tube on the open end. The time between photographs is 
given by Ai. 

The flames in the figure 1.4 take on a variety of shapes. 29 The hook-shaped flame, as seen in the 
third tube from the top, is due to buoyancy effects. A tulip flame can appear in this type of tube 
as well, as seen in the second and twelveth tubes, but it is not always present. The existence of the 
tulip in this configuration was also shown numerically. 16 The flame can also take on a ‘mushroom’ 
or parabolic shape. 

Lee and Tsai 33 and Hackert et al. 34 performed numerical studies on the effects of the wall 
boundary condition on the flame shapes, with the objective of determining whether a tulip or 
a mushroom flame is formed. For insulated tubes, both shapes appeared in tubes with a large 
diameter. As the diameter decreases, only the tulip flame shape develops. 33 The tulip flames also 
propagate at a faster speed then the mushroom shaped flames. 33 For an isothermal wall, both flame 
shapes exist, but as the diameter decreases, only the mushroom shape occurs. 33,34 For this boundary 
condition, the mushroom-shaped flame propagates the fastest. 33 If the walls are convectively cooled, 
both flame shapes appear. 34 

1.2.3 Closed - Open 

Figure 1.5 shows the channel configuration for the closed-open geometry. Here, ignition occurs 
at the closed end, and then the flame forms and propagates toward the open end, where the gas 
escapes the channel. Since the gas is allowed to escape, there is a flow formed ahead of the flame. 
This flow is caused by the expansion of the gas as it burns and pushes the flow ahead of the flame, 
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similar to a piston. Due to the motion ahead of the flame, the flame propagates at a greater speed 
then it did for the previous two channel configurations. 9 



Figure 1.5: Closed-Open Channel Configuration. 


Little experimental work 7,9,27 has been done on the closed-open configuration and no references 
to numerical work were found. Figures 1.6, 1.7, and 1.8 show schlieren photographs for flames in 
fuel lean, stoichiometric, and fuel rich conditions. 7 On the top of figure 1.6, the experimental tube is 
seen with its four windows. The schlieren photographs are shown below the window from which they 
were taken. The times on the right between the photographs represent the time intervals between 
photographs. The flame was ignited with a spark. 

After the flame induced by the spark has reached the side wall, the flames in all three mixtures 
keep their rounded shapes, no tulips forms. Guenoche reported that tulips have formed in closed- 
open channel configurations. 9 As the flames continued to propagate down the tube, they underwent 
a transition from laminar to turbulent. When this happens, the flame speed greatly increases due the 
increased burning area of the flame. The transition occurs sooner the closer that the mixture is to 
stoichiometric conditions. 7 Figure 1.8, a rich fuel mixture, shows that the flame oscillates, and even 
reverses direction for a time before resuming movement toward the open end of the tube. Mallard 
and Le Chatelier called this the ‘movement saccade’ or ‘jerky motion’. 27 
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Figure 1 .7: Flame propagation in a tube closed at the ignition end and open at the other, stoichiometric, from Schmidt et al. 




aa 


........ i ;; S 

'i 

; . : Ji H : ^ ^ : ; ; : iff 

• = . : : \ . , : 
-S- l.s' ; -='k .: •: •• S ■:-.' . 7 




4; iS ; ' : - ■ ■ N4 - ; :h ■ . ' : ii;i ! :; 

■. -■ ■■- -Si.; ■} r -77.7 

, ' v y f - : : '7 rr 

•"• ■•’ a' M:. Wl: '•»:•< ' i.- v S:!'- • 

s ■ , | . ;: • ;;■: : W : § 

■ ■ ■•; ■ ■ .' 4 - >j 7 ' '. : .■ 7 Ji 

■ ■■ . : "■ v 

: : u 7 . 


■ . ; ^'>73 ':- 

, : | ; :|; si 1 1: > iiiiii 

• i : • ■ ••" ■■■ V' { - : ' - 


" :! - H4.4 

: 774 77 7 , 

; : :S47777f7S777 ’ ; : -r : '-:f :: 

• i. : .••■.. 7b4<< ;'1:7..7:.:77: S.: : ^ : ^^xv^/iSSj : ■ 

: '•■ : • ' :• . ' ' • . 

-;is ?=r.aa:i . 

^ :: ; «; 


• Wi : iiiliiiiiiis: 

. :, 77 ' 

«E 1:[' -S’:: 




Figure 1 .8: Flame propagation in a tube closed at the ignition end and open at the other, fuel rich, from Schmidt et al. 
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Chapter 2 

Governing Equations and 
Numerical Procedure 


This chapter details the governing equations, chemistry model, integration technique, and bound- 
ary conditions. 

2.1 Navier-Stokes Equations 


The unsteady, two-dimensional, compressible, chemically reacting, Navier-Stokes equations in 
conservation form were solved. These equations model the conservation of total mass density, species 
mass in terms of number densities, momentum, and energy. In differential form, they may be 
written : 35 


• Conservation of Total Mass 

dp dpu dpv 

dt dx dy 

• Conservation of x-Momentum 

djpu) d(puu) djptw) _ dp dr xx dr zy 

dt dx dy dx dx dy 

• Conservation of y-Momentum 

djpv) djpvu ) djfiVV) _ dp dTyx cfryy 

dt dx dy dy dx dy 

• Conservation of Total Energy 

d(pet) djpetu ) djpetv) _ djpu) d(pv) dq x dq y 

dt dx dy dx dy dx dy 


(2.1) 


( 2 . 2 ) 


(2.3) 


(2.4) 


, d(UT XZ + VT X y) d(UT yX + VTyy) 

dx dy 
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• Conservation of Species Mass Density 

dm , d( ni u) , d(n iV ) d(mU Xi ) d(mU yi ) . . 

~dt + + ~ dx fy~ +UJi (2 5) 

The total internal energy is the sum of the individual species sensible internal energies plus zero 
point energies, as represented by the heats of formation, and the kinetic energy. The total internal 
energy density can be written as 

I., 

P e t = ( e,i + ^/i) + 2 P ( u2 + v2 ) (2-6) 


where ns is the number of species. The viscous stress terms, r xx , r xy = r yx , and r yiM are written as 



x ( du dv\ du 

(2.7) 

T xx 

- X \di + d^) + 2f *di 


(du dv\ 

(2.8) 

T xy 

II 

ST 

+ 

st 


x ( du dv\ dv 

(2.9) 

Tyy 

- x \di + ai) + **5 

• It . - . O 


Stokes Hypothesis is used to relate the bulk viscosity to the dynamic viscosity, A = — |/i. The heat 
flux terms, q x and q v , contain contributions due to conduction and heat flux due to the diffusion of 
the species, 


Qx = ~k^ + Y^ p iU Xi (h ti + Ah fi ) (2.10) 

dT 

Qv = + 5Z Pi U Vi +&h fi ) (2.11) 


The diffusion velocities, U x . and U y{ are calculated using a modified form of Fick’s Law, 36 


tt i-n ^Ci dCj 

e,U,. = -,d,— + 

>=1 

(2.12) 

TT r-w d&i dCj 

Pi U Vx - —pDi — + pi ^ Dj — 

(2.13) 

To complete the system of equations, the individual gases are assumed to be thermally perfect, and 

the perfect gas equation of state is used, 


P = pR*T 

(2.14) 

e »v — e s< CO 

(2.15) 


The exact formulation of the species sensible internal energy will be given in the next section. The 
species densities are constrained such that the sum of the species densities must equal the total 
density. 


13 



2.2 Chemistry Model 


A single-step chemistry model consisting of two species was implemented. The model was de- 
veloped by Khokhlov et al. 37 to simulate a stoichiometric acetylene-air reaction at 100 torr (100 torr 
= 1/7.6 atm). 

C 2 H 2 + \ {0 2 + 3.76 N 2 ) -► 2 C0 2 + H 2 0 + 9.4iV 2 (2.16) 

The model consists of a unimolecular reaction that consumes a fuel, F, representing the acetylene 
and air mixture, producing a product, P, representing the carbon dioxide, water, and nitrogen. 


F-+ P 


(2.17) 


The rate equation is written in terms of the mass fraction of fuel, Y. It is first order with the rate 
constant given by an Arrhenius Law, 


dY 

dt 


= —pYkf 


(2.18) 


k f = A T e ( - T ' /T) (2.19) 

The energy release due to the combustion is a constant given by 

q = 35 R s T r (2.20) 

This represents an effective zero-point energy or heat of formation of the fuel species, F. The heat 
of formation of the product species, P, is zero. The use of a simplified one-step chemistry model, 
as opposed to a multi-reaction, multi-species model, is done to reduce the chemical integration time 
while still maintaining an accurate representation of the chemistry in the problem. A multi-species, 
multi- reaction model can lead to significant increases in the computational time. In some situations, 
the chemistry integration time can be of the same order of magnitude or greater then the fluid 
dynamic integration time. 38 

Table 2.1 gives the values of the chemical and gas properties. The molecular weights of the fuel 
and product are equal and constant. The value of 29.0 was chosen as the average of the molecule 
weights of the acetylene and air mixture and of the carbon-dioxide, water and nitrogen mixture. 
The specific heats, Cp and Ct,, of both species are equal and constant. The actual values of the 
specific heats increase as the temperature increases. Since one of the goals of this research is to 
investigate the basic flow physics, the use of a constant specific heat assists in the formation of 
constant flow parameters, such as Prandtl number and Lewis number. Also, the use of a constant 
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Table 2.1: Chemical and Fluid Properties 


Property 

Value 

M 

29.0 (g/mol) 

R. 

2.8670 x 10 6 (erg/g-K) 

7 

1.250 

°v 

1.4335 x 10 7 (erg/g-K) 

Cv 

1.1468 x 10 7 (erg/g-K) 

A r 

1.0 x 10 12 (cm 3 /g-s) 

T r 

293 (K) 

T a 

29.37V (K) 

Hr 

2.439 x 10" 8 (g/cm-s) 

k r 

2.439 x 10" 8 (g/cm-s) 

D r 

2.439 x 10 -8 (g/cm-s) 


specific heat reduces the cost of the computation by simplifying the calculation of the temperature. 
If a function of temperature was implemented, the temperature calculation would be an iterative 
procedure, increasing the computational time. These parameters were chosen by Khokhlov et al. 37 
to fit the measured values of the laminar flame speed, laminar flame thickness, Chapman- Jouget 
detonation velocity, and the detonation thickness of an acetylene-air mixture for a pressure range of 
0.1 to 1.0 atmospheres. 

With the specific heats constant, the sensible internal energy of both species becomes c^T. It is 
now possible to define the equations for the species interned energies and enthalpies. Both contain 
the sensible part and the heat of formation. 


ey = CvT + q 

(2.21) 

ep ~ CvT 

(2.22) 

h Y = CpT + q 

(2.23) 

hp = CyT 

(2.24) 

Taking the above equations etnd substituting them into equation (2.6) 
total internal energy as 

gives the equation for the 

pe t = pevT + (u 2 + v 2 ) + pYq 

(2.25) 
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2.3 Transport Properties 


2.3.1 Dynamic Viscosity and Thermal Conductivity 


The dynamic viscosity, /i, and the thermal conductivity, Jk, are calculated from power-law func- 
tions of temperature, 37 



(2.26) 

(2.27) 


where n = 0.7. These forms of the equations give a Prandtl Number of 


Pr = ^2 = £ (2.28) 

k k r 

This allows parametric studies to be done by changing the reference coefficients, p r and k r , to get 
the desired Prandtl Number. 


2 . 3.2 Mass Diffusion Coefficient 


Since in this model the gas contains two species, the binary diffusion coefficients for the two 
species are equal, 39 Dy = Dp. These coefficients are calculated as power-law functions of tempera- 
ture and density, 37 



P 



(2.29) 


With the binary diffusion coefficient in this form the Lewis-Semenov number can be written as 


Le = = — 

k k r 


(2.30) 


This definition of the Lewis number is the standard way it is defined in aerodynamic texts. Some 
texts on combustion define the Lewis Number as the reciprocal, so that Le would be equal to k r /D r . 
For a binary gas, Fick’s Law as given by equation (2.12) reduces to 


P,v„ = -pD a £ 


PiU Vi = -pD 


dcj 

dy 


(2.31) 

(2.32) 


Using this form of Fick’s Law, the definitions of the enthalpies as given in equation (2.21), and the 
fact that the sum of the mass fluxes must be zero, (pyUy = — ppUp ), the heat flux terms from 
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equation (2.10) can be written as 


(2.33) 

(2.34) 


,&T _ dY 

= -* 5 ? - ' ,D,! ^ 

.«r ar 

* = -* a? “ ^ to 


2.4 Numerical Procedure - LCPFCT 


The convective terms of the Navier-Stokes equations are integrated in time using a parallel 
version of the Flux Corrected Transport model, LCPFCT. 40 The viscous, conductive, and mass 
diffusion terms can be treated as either source terms or can be added in separately after the convective 
terms are integrated. The pressure terms are always included as source terms for the LCPFCT 
routine. The species production term, a)*, is integrated separately. 

The LCPFCT routine is a cell centered, finite volume convection routine which is second-order 
in time and fourth-order in space. The evaluation of the source terms, the viscous, conduction, and 
mass diffusion fluxes, are second order, which reduces the entire scheme to second order in space. 
LCPFCT is a nonlinear, monotone method designed to calculate sharp gradients with minimal 
numerical diffusion. The flux quantity, for example p, is convected to an intermediate value in a 
predictor step. Diffusion is added to guarantee that the new density is monotonic. This value is 
then corrected by an anti-diffusion term to remove the added diffusion and advance the density to 
the new time. The amount of anti-diffusion that is added depends on the local values of the flux 
term making the method nonlinear. The amount of anti-diffusion added also ensures that no new 
minima or maxima are generated, or adds to those that already exist. 

The LCPFCT routine solves the general, one-dimensional, continuity equation as given by : 


_ 1 d(r° 1 jm) 1 d(r a ~ l Bi) dB 2 _ 

dt r a-l Q r r o-l Q r 2 dr 3 


(2.35) 


The density flux terms in the Navier-Stokes equations, p, pu, pu, and pe t , are represented by ip. The 
a represents the coordinate system used, 1 for Cartesian, 2 for cylindrical, and 3 for spherical, with 
r being a generalized coordinate direction. The terms B \ , B 2 , and Bz represent source terms, such 
as viscous stress terms or heat flux terms. The C 2 is a constant. 

LCPFCT first convects the the above general continuity equation to get a temporary term ip* 


A Wi = A ° x xP°i - Au i+J + Atip°_ h (2.36) 

where A is the cell volume per unit depth, A is the cell interface area per unit depth, and Au is 
the difference between the flow velocity and the grid motion for a moving grid. For the stationary 
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Cartesian grid that was used for this study, A = 1, A< = - r 4 _j, and Au is equal to the fluid 

velocity. Since LCPFCT is a cell centered algorithm, the subscript i refers to the cell centers. The 
subscripts i ± 1/2 refer to points on the cell interfaces. These cell interface values are calculated as 
average values of the cell centered values that they are located between. The source terms are then 
added in to calculate a transported density, rj)J. 

A WT = + \±tA i+i (B u+1 + B u ) - (B lti + *!,*_,) 

+ -A<C2,i + A-i) (i?2,»+i - f?2,«-i) (2.37) 

The diffusion is now added to obtain the transported-diffiised density, tpi. 

A itii = Ktf + K+\ Ai+i (V'f+i - O ~ "i-jAi .1 (V>? - Vi-x) ( 2 . 38 ) 

where the diffusion term, i/, is given by 

v i+l = l + hu (2- 39 ) 

€ i+i =.A 1+4 4« j+ j T (- + — j (2.40) 

The anti-diffusion flux, , is calculated from the transported density, , from equation (2.37). 
This density is used so that there is no residual diffusion for the case when the fluid velocity is equal 
to the grid velocity (Au i+ ^ = 0). This would be the case if the transported-diffusion density, xp, 
was used instead. 

/“+! = /*i+iA i+i (rpf +1 - tf) (2.41) 

where the anti-diffusion term, jx, is given by 

m l~ Hi (242) 

The diffusion term, v, and the anti-diffusion term, (i, were chosen by Boris and Book 41 to reduce 
the phase errors for the convection to be fourth order accurate. The anti-diffusion flux must now be 
corrected to ensure that no new minima or maxima are created. This is accomplished by applying 
a flux limiter to obtain the corrected anti-diffusion flux, f c , . 

fi+\ ~ mln [i/i+jl> Aj + 1 2 — f/’i+l) ,S i+ ^A” (t/'i — V’i-l)] 

f- +i = S i+i max{0,/; +i } (2.43) 
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where S i+i is given by 


= 5 x sgn (^+1 - rp^ (2.44) 

In theory S should be equal to one. In practice the value used is slightly less then one, around 
0.999 to 0.997, to provide residual diffusion to be retained by the solution. Then, from the corrected 
anti-diffusion flux, the new density, xj )\ \ is calculated by 

(2 ' 45) 

To obtain second order accuracy in time, the density flux terms from the Navier-Stokes equation, 
P > pu, pv, pe t , and n, are convected to the half time-step, t? = t° + t. These new half-step values 
are used to evaluate the density fluxes for the entire time step, t n = t° + At. 

To use the one-dimensional LCPFCT on the two-dimensional Navier-Stokes equations, the equa- 
tions are split into two sets of equations, one for the x direction and the other for the y direction. 
Each direction is then integrated separately. To avoid biasing one direction over the other, the inte- 
gration for each time step is switched. For one time step the x direction is done first and then the 
y, for the next time step the y is integrated first then the x. To perform this direction splitting, the 
time step must be small enough so that the cell values do not change rapidly over the time step. 40 

2.5 Time Step Determination 


The time step for explicit methods, such as LCPFCT, axe restricted due to stability requirements. 
These can be derived as in Hoffman 42 for the linear stability of parabolic equations. The time step is 
restricted to be the minimum value of the convective, viscous, conductive, and diffusive time steps. 


At (CFL) min (At conu , A£ V j SC , A t con< j, A £<ft//) 


(2.46) 


where CFL represents the Courant-Friedrichs and Lewy number. For each new time step, the above 
equation is evaluated over the entire grid. For positivity of the LCPFCT routine, the CFL number 
is limited to be less then one-half. 40 It was determined that for the flame calculations, the CFL 
number needed to be less then 0.25 for stability. A value of 0.2 was used for all calculations. 

The convective time step is determined by 


A^conv — TTIXTI 

The viscous time step is determined by 


VKil + c*,/ Kil + <W 


— TTiin 


A yjA 

K 2l /jj J 


(2.47) 


(2.48) 
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where v is the kinematic viscosity coefficient. The thermal conduction time step is determined by 


Atcond — 7TIZTI 



where a is the thermal diffusivity. The time step due to the mass diffusion is given by 


(2.49) 



(2.50) 


where D is the binary mass diffusion coefficient. The limiting time scale was determined to be the 
convective time-step in the burned region. 


2.6 Chemistry Integration 


After LCPFCT has solved the convective part of the Navier-Stokes equations the chemistry 
part is integrated. The chemical source term, a), appears only in the species continuity equation 
(2.5). The other conserved quantities, p, pu , pv, and pet, are held constant during the chemistry 
integration. The chemical reactions do not change the total internal energy, pe t , as given by equation 
(2.25) but only convert the chemical energy, pY q, into sensible energy, pc^T . 

From the chemical rate equation, (2.18), the rate constant defined in equation (2.19) is evaluated 
at the current temperature. The resulting ordinary differential equation is then integrated in time 
to give the new species fuel mass fraction as 

Y” = Y°exp{—pkfAt) (2.51) 


With the new fuel mass fraction found from the above equation, the product mass fraction is found 
from P = 1 - Y since the sum of the mass fractions must equal one. The new number densities are 
calculated from the mass fraction of fuel as 


pYN a 

nY = - at 

p( 1 - Y)N a 
np= M 


(2.52) 

(2.53) 


The temperature and pressure are updated due to the reaction. EYom the equation for the 
definition of the total interned energy (2.6) the temperature can be found as 


pet ~ Ip (u 2 + 1 > 2 ) - pY^q 

PCv 


(2.54) 


The pressure can be found from the equation of state (2.14). 
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2.7 Boundary Conditions 


The Navier-Stokes equations control the fluid dynamics of the system. The boundary conditions 
are what defines the system. These conditions tell the equations where there are walls, inflows, 
outflows, etc. The boundary conditions control the solution to the problem. 

The boundary conditions for subsonic inflows, subsonic outflows and walls were calculated using 
the characteristic method of Poinsot and Lele 43 developed for a three-dimensional calorically perfect 
gas. Their method for boundary conditions of the Navier-Stokes equations is an extension on the 
Euler method of Thompson . 44,45 This method uses the Navier-Stokes equations to solve for the 
boundary variables that are not specified. It does this by calculating the characteristic waves that 
are leaving and entering the system and using these waves to calculate the variables that are not 
known. 

Using the characteristic analysis in Thompson , 44 the Navier Stokes equations are written for a 
two-dimensional system as 43 : 


dp . d(pv) 

m +d ' + -dT 


= 0 


(2.55) 


djpv Q 

dt 


+ ud i + pdz + 


d{puv) 

dy 


dT X y 

dx “aiT 


(2.56) 


djpv) 

dt 


4- vd\ 4- pdi + 


d(pvv) 

dy 


dp djj/X dTyy 

dy dx dy 


(2.57) 


^ {u 2 + v 2 ) dr + + pud 3 + pvd 4 


^ djpetv) _ djpv) dqx dqy 

dy dy dx dy 


d(uT xx + VT xy ) d(ur yx + VTyy) 

T /-v 


(2.58) 


dx dy 

It is assumed that the x-direction is normal to the boundaries. The dj terms contain the values for 
the derivatives in the x direction. They are given by : 44 


di 

^3 

d 4 


c 2 


*2 + 2 (*1 + *4) 


-(*l + * 4 ) 

¥3 


(2.59) 


21 


(2.60) 

(2.61) 

(2.62) 



inflow 


outflow 


^4 

*3 

*2 




flow direction 


*4 

*3 

*2 


Figure 2.1: Characteristic Wave Configuration 




x-axis 


The terms are the amplitudes of the characteristic waves. These waves are shown in Figure 2.1 
for an the inflow and outflow configuration. The waves $2 , ^3, and ^4 move in the positive x 
direction. The wave $1 moves in the negative x direction. These terms are given by : 

*1 

*2 

*3 


The waves propagating from the inside of the domain out of the boundaries are calculated 
using one sided finite differencing. To calculate the waves moving into the system from the outside, 
Poinsot and Lele 43 used a one-dimensional analysis to estimate their values. In terms of the primitive 
variables, p, u, v, p, and T, these equations are written as 44 : 



(2.63) 


(2.64) 

( dv\ 

u \fo) 

(2.65) 

, . (dp 6u\ 

{ „ +c) i- +pc -\ 

(2.66) 


dp 

£ +d '=° 

du 

m +d, = 0 
£ + 4 = 0 

dp 

£+4=0 


dT T_ 
dt p<? 


-*2 + 5(7 - 1)(*1 + *4) 


= 0 


(2.67) 

( 2 . 68 ) 

(2.69) 

(2.70) 

(2.71) 


Prom these equation, relations between the incoming and outgoing waves can be established based 
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on the set values at the boundaries. 

Saint-Martin-Tillet 46 implemented these boundary conditions for nonreacting flow in a two- 
dimensional channel. He modeled the boundary conditions by considering the one-dimensional 
Euler equations alone as given by equations (2.67 - 2.71). This worked for the wall boundary and 
for the inflow if it consisted of parallel flow, v = 0. It also worked well for the outflow in an inviscid 
region of the flow, such as near the centerline of a channel. Problems arose for inflows that were 
not parallel, such as flow at a stagnation point, and for the outflow of a viscous region, such as 
a boundary layer. These problems were corrected by expanding the boundary conditions to two 
dimensions and solving equations (2.55 - 2.58), including all the viscous terms. 

2.7.1 Subsonic Inflow Boundary Condition 

For an inflow boundary, it is common to define the velocities and a thermodynamic variable, 
such as temperature. For a reacting case, the inflow mass fractions are also specified. With the 
primitive variables u, v, and T specified the boundary equations (2.56), (2.57), and (2.58) do not 
need to be solved, only the equation of the global mass density, (2.55) needs to be solved. 

For an inflow, the wave *i is leaving the domain and the other waves are entering, therefore *i 
is the only wave that can be calculated. To obtain values for the waves that are entering the domain 
through the boundary, the one-dimensional equations (2.67 - 2.71) are used to define the incoming 
waves in terms of the outgoing wave. Imposing the velocities (u and v are set) gives d$ = 0 and 
d± = 0. Using equations (2.61) and (2.62), these conditions define the the wave amplitudes *3 and 
$4 as 


*3 = 0 (2.72) 

* 4 = *1 (2.73) 

From equation (2.71), imposing the temperature and making use of the above equations gives the 
expression for *2 as 


* 2 = ( 7 - l)*i (2.74) 

These values are used to calculate d\ which is used in the calculation of the new inflow density. Once 
the new global density is determined the pressure is obtained from the equation of state. The inflow 
number densities are calculated from the specified inflow mass fraction and the new global density. 
The total energy density is calculated from equation (2.25). 
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2.7.2 Wall Boundary Conditions 

For the case of a solid wall, the no-slip condition applies. This means that the two components 
of velocity are zero at the wall, u = v = 0. Using u = 0 in equations (2.63) and (2.65) gives the 
waves ^2 and $3 as 

*2 = *3 = 0 (2.75) 

From the one-dimensional equation for u, (2.68), <*3 = 0 which provides a relationship from the 
remaining two waves as 

! = *4 (2.76) 

From Figure 2.1, if the inflow is replaced by a wall, then the wave $1 is calculated from the compu- 
tational system. If the outflow is replaced by a wall then the wave ^4 is calculated. 

For the wall temperature, there are three possible numerical boundary conditions. The first is 
the isothermal wall where the wall temperature is held constant. The second is the adiabatic wall 
where the heat flux to the wall is zero. The third possible numerical boundary condition for the 
temperature at the wall is the convectively cooled or heated wall. The wall temperature is calculated 
based on an energy balance between the wall heat flux and the convective heat flux. If the wall has 
a non-zero thickness, then modeling of the heat conduction through the wall is required. This third 
type will not be discussed here. 

For the mass fraction at the wall there are three possibilities. The first is an equilibrium wall 
where the mass fractions are calculated based on their chemical equilibrium value at the wall tem- 
perature. The second is a non-catalytic wall, where the mass flux to the wall, as given by Fick’s’ 
Law, is zero. For the binary gas system this condition reduces to a zero mass fraction gradient at 
the wall. The third condition, lying in between the first two conditions, is for a catalytic wall. For 
this boundary condition the species mass fractions are calculated based on a combination of the wall 
material properties and mass diffusion flux from the fluid. This third type will not be discussed 
here. For an excellent explanation of catalytic wall boundary conditions see Grumet. 47 

Using the first two temperature and mass fractions boundary conditions, four sets of wall condi- 
tions can be formed. 

• For an isothermal, non-catalytic wall, the temperature is specified and the species mass frac- 
tions are found by setting the mass fraction gradients at the wall to zero. 

• For an isothermal equilibrium wall, the temperature is specified. The species mass fractions 
are calculated based on the wall temperature. 
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• For an adiabatic, non-catalytic wall, the mass fraction is found by setting the gradient at the 
wall to zero. The temperature is found by setting the temperature gradient at the wall to zero. 
This has the effect of making the heat flux to the wall zero. 

• For an adiabatic, equilibrium wall, the wall conditions need to be calculated based on the 
condition of zero heat flux from equation (2.33). Assuming the wall is normal to the x direction 
this gives 

<277) 

As a first estimate, the the old values of the wall mass fraction sire used to calculate the 
species gradient. The above equation can then be used to iterate on the new wall temper a- 
ture. An iteration method is required since the species mass fraction is a function of the wall 
temperature. 

With the wall temperature, mass fraction, and velocities known, only equation (2.55) for the 
density at the wall needs to be solved. With the new wall density, the equation of state is used 
to calculate the wall pressure. The definition of mass fraction is used to calculate the new number 
densities. The total energy density at the wall is calculated from equation (2.25). 

2.7.3 Subsonic Outflow Boundary Conditions 

For the subsonic outflow, pressure waves can both enter and leave the computational domain 
through the boundary. The wave, that is entering the domain carries information into the 
domain from the downstream conditions. The other waves, $3, and $ 4 , are leaving the system 
and can be calculated from the interior grid points. The wave entering the system needs to be 
estimated. Pomsot and Lele did this by specifying a constant pressure at infinity and calculating 
the entering wave as 43 


*x = K(p- Poc ) 


(2.78) 


The constant, K, taken from Rudy and Strikwerde, 48 is given by 


K = 


ac( 1 — M ) 


(2.79) 


where M is the maximum Mach number in the flow, L is a characteristic length such as channel 
height or plate length, and a is a constant equal to 0.25. 43 

For a viscous problem, there are additional viscous boundary conditions that axe required. 43 
These viscous boundary conditions for a two-dimensional flow axe that the derivative normal to the 
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boundary of the tangential shear stress, r zy and the normal heat flux, g z , are zero. These conditions 
are implemented directly into equations (2.55 - 2.58) when they are solved to get the outflow fluxes. 
The boundary condition for the species are obtained from simple extrapolation of the mass fractions 
for the interior. 

Calculating the outflow boundary conditions using the above methods works for single and mul- 
tispecies nonreacting flows. When the chemistry was turned on and the flow allowed to react the 
outflow boundary conditions became unstable when the flow was treated as two-dimensional. For 
the one-dimensional reacting case the boundary conditions described above were stable. To alleviate 
the unstable problem the energy density and mass density were extrapolated at the outflow bound- 
ary. The velocities were calculated using equations (2.56) and (2.57) so that the characteristic waves 
could be accounted for. The species mass fractions were also extrapolated. One possible reason for 
the instability is that the boundary condition model was developed for a non-reacting, calorically 
perfect gas. 

2.8 Initial Conditions 

The initial channel configuration is shown in figure 2.2. The channel is 8.2 cm in length and 
0.127 cm in width. The left and bottom sides are solid walls. The top of the domain is a symmetry 
condition and the right side is an outflow. This simulates a rectangular channel that is closed on one 
end and open at the other. The symmetry condition is used here since the gravitational terms have 
been neglected. If the gravitational terms were included, then the orientation of the channel would 
be important. If the tube was horizontal, buoyancy effect would distort the flame and a symmetry 
condition could not be used. If the channel was vertical, then a symmetry could be used. 

Starting at the closed end, the first 0.3 cm is set to the adiabatically burned conditions. These 
conditions come from laminar flame theory that will be discussed in section 4.2. From 0.3 cm to 
the end of the channel, the gas is at the unburned conditions. The values of the properties in the 
burned and unburned region are given in Table 2.2. 
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Table 2.2: Initial Conditions 


Value 

Unburned 

Burned 

T (K) 

293 

2344 

p (dynes) 

133322 

133322 

P (g/cm 3 ) 

1.587 x 10~ 4 

1.984 x 10 -5 

Y 

1 

0 

u (cm/s) 

0 

0 

v (cm/s) 

0 

0 



h 


8.2 cm 


Figure 2.2: Two-Dimensional Channel Configuration 
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Chapter 3 

Code Validation 


The present chapter consists of three sections, each with various subsections. The first section 
describes previous uses of the parallel code. The second section contains results for a non-reacting 
flow. The third section incorporates the chemistry with the fluid dynamics. 

3.1 Past Code Work 

CMRFAST2D was used successfully in previous works. J. Weber 38 parallelized an Euler code and 
added chemical reaction to study the two dimensional detonation wave problem using a hydrogen-air 
chemistry. Y. Weber 49 added viscosity and thermal conductivity to the Euler code and calculated 
the interaction of a reflected shock with a boundary layer in a shock tube. Piana 50 added the effects 
of mass diffusion to create a fully reactive Navier-Stokes code. Nguyen et al. 51 used a non-reacting 
version to study the interaction of a vortex with a boundary layer in a channel containing ribs. 

3.2 Non-Reaction Code Validation 


The non-reacting parts of the code were tested to ensure that they are working properly. These 
parts consist of the viscous, thermal conduction, mass diffusion, and Euler convection routines. 

3.2.1 Gaussian Profile - Mass Diffusion 

The convection, viscosity, and thermal conductivity are turned off. The pressure and tempera- 
ture are held constant at 100 torr and 293 K respectively, which fixes the mass diffusion coefficient 
and global density. This reduces the governing equations to just the two species continuity equations. 
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Table 3.1: Gaussian Profile Parameters 


Parameter 

Value 

P 

1.723 x 10" & (g/cm 3 ) 

D 

0.402 (cm 2 /s) 

U 

0.15 (s) 

n i 

7.876 x 10 17 (1/cm 3 ) 


Making use of Fick’s Law and reducing to one dimension, equation (2.5) becomes : 

drii d 2 n< 

dt ax 2 


(3.1) 


Initially, species 1 is deposited along the x axis from -0.1 to 0.1 at a value equal to the global 
density, and is set to zero everywhere else. This species density spreads out due to the diffusion and 
becomes Gaussian over time. The theoretical solution for the density profiles in time and space for 
the Gaussian diffusion problem is given by: 52 


n(x,t) = ni 




(3-2) 


In this equation, t\ and ni refer to a reference state that occurs after the initial start of the diffusion 
process. For this problem, the reference state was chosen at 0.15 s. The number density for the 
reference state at x = 0 is n\ = n(0, ti). Table 3.1 gives the values of the parameters used for 
this calculation. Figure 3.1 shows the numerically generated Gaussian profiles versus the theoretical 
values at 0.3 s and 0.5 s for the first species. The symbols represent numerical values and the lines 
represent theoretical values. The cell size is constant at 0.02 cm. The code and theory agree well. 


3.2.2 Heat Flow Problem - Thermal Conductivity 


For this problem the viscosity, convection, and mass diffusion are turned off. The gas used was 
air with a constant density of 1.548 x 10 -4 g/cm 3 which is based on a pressure of 100 torr and a 
temperature of 300 K. The thermal conductivity was also set constant, at 2650 erg/cm-s-K, which 
was calculated from Sutherland’s Law 53 at 300 K. With these assumptions, the energy equation (2.4) 
in one dimension reduces to : 

,33) 

dt~\pc,)dx‘ ' 

This is the classical form of the heat conduction equation. 

Initially, the temperature of the domain was constant at 300 K. Then the side at x = 0 was 
increased impulsively to 600 K and held at that temperature. The classical solution of the conduction 
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x (cm) 

Figure 3.1: Gaussian Profiles for times 0.3 s and 0.5 s. 
equation from Hildebrand 54 for these conditions reduces to : 

T(x,t) = [7\ + (T 2 + TO J (3.4) 

nn L 

n=l 

where T\ = 600 K, Ti = 300 K ,1 = 1.0 cm, and fc is given by : 



Figure 3.2 shows the result of the integration for a time of 12.2 ms. The grid cell size is 0.016 cm. 
The solid line is the theoretical value and the symbols are the numerical values. The summation in 
the theoretical equation was taken to n = 600. The plot shows an excellent agreement between the 
theoretical and numerical temperatures. 

3.2.3 Stokes’ First Problem - Viscosity 

For this case, the convection, thermal diffusion, and mass diffusion are turned off. The density, 
pressure, and temperature are all held constant. The viscosity coefficient is also held constant at 
1.846 x 10“ 4 g/cm-s, which is the value at 300 K from Sutherland’s Law. 53 The density is 1.548 x 10 -4 
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Figure 3.2: Temperature Profile for the One-Dimensional Thermal Conduction Equation. 


g/cm 3 . In one dimension, the momentum equation reduces to 


du _ 4 d 2 u 
dt 3 dy 2 


(3.5) 


Initially the air in the domain was at rest. Then at y = 0, the velocity w’as impulsively increased 
to 600 cm/s and held at this speed. The classical solution to this problem is given by Schlichting 55 
as 


u _ 2 f° 

U.-y/iL 


exp(—T} 2 )dr) 


where U 0 = 600 cm/s, and rj is the normalized length scale given by 

V 


r) = 


2 y/ut 


(3.6) 


(3.7) 


Figure 3.3 shows the normalized calculated values as compared to the classical solution. The 
agreement is excellent. The grid resolution for this problem is 0.006 cm. An interesting result of this 
problem is that the boundary layer thickness is proportional to the square root of the product of the 
kinematic viscosity and time, <5 oc \fvi. This result, and the results for varying edge flow velocity in 
Appendix A, will be useful in the discussion of the formations of the boundary layer in the channel. 
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Figure 3.3: Normalized velocity profile for Stokes’ first problem. 

3.2.4 One Dimensional Shock Wave - Convection 


For this case, the viscosity, thermal conduction, and mass diffusion are turned off. This reduces 
the Navier-Stokes equations to the Euler equations. A one dimensional shock tube problem for air 
was set up with an initial diaphragm pressure ratio, Pa/P\, of 120. The initial temperature was 
uniform at 300 K. The theoretical solution to this problem can be found in Anderson. 56 From the 
diaphragm pressure ratio, the pressure ratio across the shock, pz/pu can be calculated from 

P4_P2[, (7-l)((P2/Pi)-l) 1'^ 

Pi Pi \ i/ 2 7[27 + (7 + l)((P2/Pi) - 1)] J l J ° 

where the conditions in region 2 are behind the wave and the conditions in region 1 are ahead of the 
wave. Once the pressure ratio is calculated, the density and temperature ratios across the shock are 
calculated from 



Tl — El/£l 

T\ Pi Pi 


(3.9) 

(3.10) 
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Table 3.2: One-Dimension Shock Properties 


Value 

Code 

Theory 

P 2 /P 1 

6.75 

6.75 

P 2 M 

3.26 

3.25 

t 2 /ti 

2.07 

2.07 

JfE 

585.1 m/s 

584.0 m/s 


The piston velocity, u p , is the velocity that the shock wave induces into the flow behind it. It is 
calculated from 



(3-11) 


Figure 3.4 shows the shock pressure ratio, density ratio and piston velocity. Table 3.2 gives the 
calculated shock properties versus the theoretical calculated values. The grid resolution for the 
numerical solution is 0.03 cm. As can be seen, all the code values compare to well with 1% of the 
theoretical values. 


3.3 Reacting Code Validation 


For these tests the chemistry routine was active in the integration. The chemistry routine was 
tested for a constant volume combustion and a one-dimensional detonation wave. The results for a 
one-dimensional laminar flame will be presented in the next chapter. 

3.3.1 Constant Volume Combustion - Chemistry Alone 

For this case, the convection is turned off along with all the transport properties. The mass 
density, p, and the total energy density, pe f , are constants. The energy equation then reduces to 

CyTf + qYf = CvTi + qYi (3.12) 

The initial conditions for temperature and pressure were 500 K and 100 torr respectively. The initial 
mass fraction, Yi, was set to 1. Assuming that the reaction goes to completion, Yf = 0, the above 
equation gives a final temperature of 3063.75 K. By using the perfect gas law, and the fact that 
density is constant, the final pressure becomes 8.17 x 10 5 torr. 

Figure 3.5 shows plots of the temperature and pressure versus time. The code generated final 
values of temperature and pressure are 3063.75 K and 8.17 x 10 5 torr. These match the expected 
values exactly. This test shows that the chemistry integration routine is working properly. 
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Figure 3.4: Shock wave profiles for air in a shock tube with initial pressure ratio of 120. 
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Figure 3.5: Temperature and pressure for a constant volume combustion. 
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Table 3.3: One-Dimensional Detonation Wave Properties 


Values 

Theory 

Code 

% Difference 

P2/P1 

19.04 

18.30 

3.89 

pi!p\ 

1.76 

1.71 

2.84 

Ta/Ti 

10.83 

10.70 

1.20 

u p (m/s) 

808.3 

775.0 

4.12 


3.3,2 Detonation Wave - Chemistry and Convection 

This case is similar to the 1-D shock wave done above except now the chemistry integration 
routine is active. The viscous, thermal conductivity, and mass diffusion are all deactivated, leaving 
active just the chemistry and convection routines, and reducing the equations to the reacting Euler 
equations. The initial conditions were set up similarly to the shock tube problem. A diaphragm 
was placed at x = 6 cm with a pressure ratio across the diaphragm of 150 and a temperature ratio 
of 4. Burned gas, Y = 0, was placed from 0 cm to 6 cm. Unbumed fuel, Y = 1, was placed from 6 
cm to the end of the tube. The solution to the detonation wave problem can be found in numerous 
combustion textbooks. 2,5,57,58 

The governing equations in one dimension for the flow across a detonation wave are given as 


PlUl = P2U2 

(3.13) 

1 2 1 2 

Pi + 2 PlU l = P2 + 2 P2 “ 2 

(3.14) 

CpTi + -u\ +q = CpT 2 + -u 2 

(3.15) 


The equation of state also holds on both sides of the wave. 

p = pR s T (3.16) 

The condition for a Chapman-Jouget detonation is that the velocity of the burned gas, region 2, be 
at the local speed of sound. 

u 2 =a 2 = \/yR,T 2 (3.17) 

This equation (3.17) provides the closing relationship. Equations (3.13)-(3.17) provide five equations 
for the five unknowns, P2, X2, p2> Ui, and 1*2. The known values in the problem are the conditions in 
region 1, Pi , 7i , pi , the chemical energy release, <7, which is defined in equation (2.20), the specific 
gas constant, R $ , and the ratio of specific heats, 7. 

Figure 3.6 shows the computed pressure ratio, temperature ratio, and velocity of the detonation 
wave at one time. The grid resolution for this case is 0.001 cm. Table 3.3 lists the calculated 
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detonation properties, the theoretical values, and the percent difference between the two. All the 
numerically calculated values are within 5% of the theoretical values. The plot of pressure shows 
an initial spike before it relaxes to the value of 18.30. At the time shown, the spike has a pressure 
ratio of 43.9 and corresponds to the von Neumann spike as discussed in Kuo 2 and Williams. 58 The 
calculated value of the von Neumann spike, determined as discussed in Williams, 58 is found to be 37. 
The numerical value of this spike is 18% greater than the theoretical value. The magnitude of the 
spike changes in time but is always greater then the theoretical value. The longitudinal waves seen 
in the solutions are a property of “galloping” detonations. 38 For one-dimensional detonations, the 
shock front oscillates as a result of the detonation propagating into unreacted gas. The oscillations 
in shock strength caused the properties behind the shock to oscillate as well, causing the longitudinal 
waves. 
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Figure 3.6: Pressure ratio, temperature ratio, and velocity for a 1-D detonation wave. 
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Chapter 4 

One-Dimensional Flame Results 


4.1 Overview 


A series of one dimensional flame simulations were performed to evaluate the chemistry model 
with the convective and all the transport properties active. A comparison of the calculations with 
the laminar flame theory is shown. This theory is described in the next section. A grid study was 
performed to evaluate the dependence of the gridding on the flame speed and flame thickness. The 
transport properties were varied to compare the results to the general trends as provided by theory. 
The initial conditions are the same as the initial conditions described for the two-dimensional channel 
in section 2.8. For this case though, a symmetry condition is applied at the bottom boundary to 
create a one-dimensional configuration. The initial conditions axe provided in table 2.2. 

The general trends for the flow properties through a one-dimensional flame are shown in Fig- 
ure 4.1 from Kanury. 59 In the figure, the flame is moving to the right into a stagnant gas. As 
the flame propagates into the unbumed stagnant gas, the temperature increases and the density de- 
creases. The pressure slightly decreases, which is the basis for may numerical flame studies occurring 
at an assumed constant pressure. 

4.2 Laminar Flame Theory 


Laminar flame theory provides a means to estimate the flame speed, which is defined as the speed 
of the flame relative to the flow ahead of it. It also allows the estimation of the flame thickness. The 
detailed theory of Zeldovich and Frank-Kamenetskii 60-62 and Semenov 63 is described in Appendix C. 
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Prom the theory, the resulting equation for the laminar flame speed is 

-v /S)MSF 

where I is the reaction rate given as 

f T> 

1=1 udT (4.2) 

According to Kuo, 2 this formula doesn't give very accurate results, but it does predict the trends 
well. One of these trends is that the flame speed is proportional to the square root of the thermal 
conductivity, 5/ oc yfk . Another result from the theory is that for an adiabatic flame, the enthalpy 
across the flame is constant. This can be used to relate the temperature in the burned region to the 
temperature in the unburned region. 


*7 (4.3) 

CpT u + q - CyT + qY (4.4) 

Turns 64 outlines a method from Spalding 65 to calculate the laminar flame thickness, defined by 


4 = 


T b -T u 

(dT) 

V dx / max 


Ftom this method, the flame thickness can be estimated from the laminar flame speed as 

Si \PuCpJ 


(4-5) 


(4.6) 
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Since the flame speed is proportional to the square root of the thermal conductivity, this equation 
shows that the thickness is also proportional to the square root of the thermal conductivity, a y/k. 

4.3 Standard Case 

For the first case the computational cell size was set at a constant 0.001 cm, which will be shown 
later to give about 20 cells in the ID planar flame thickness. The transport properties had all the 
reference coefficients equal, as given by Table 2.1, so that Le = Pr = 1. The actual Lewis number 
for a stoichiometric acetylene and air mixture is 1.02 at 298 K and 100 torr. The actual Prandtl 
number under those same conditions is 0.75. Figure 4.2 shows the computed temperature, velocity, 
and fuel mass fraction profiles through the flame. Figure 4.3 shows the density, chemical energy 
release, and pressure through the flame. The trends show by these figures for the one-dimensional 
flame match the general trends as shown in Figure 4.1 from Kanury. 59 As the flame propagates, it 
induces flow ahead of it. The flow ahead of the flame reaches a steady condition at a velocity of 
913 cm/s, relative to fixed coordinates. The pressure curve slightly increases ahead of the flame, on 
an order of 0.05% which is consistent with the assumption that is made in laminar flame theory of 
constant pressure. 

From the temperature profile, the ratio of the temperature of the burned to unburned material 
is 8.01. From theory, the enthalpy across the flame is conserved and is given by equation (4.3). The 
right hand side of the equation is evaluated in the unburned gas, where Y = Y u = 1, and T = 
T u — 293 K. Using the definition of the chemical heat release (2.20), and noting that the reference 
temperature is equal to the unbumed gas temperature, T r = T u , the theoretical temperature ratio 
across the flame is 

^=8 (4-7) 

■L u 

This theoretical value is very close to the 8.01 obtained numerically. 

Figures 4.4 and 4.5 show the flame position and velocity in the laboratory or fixed frame of 
reference. The flame velocity, S/r, approaches 1044 cm/s relative to fixed coordinates. The laminar 
flame speed is defined as the speed of the flame relative to the flow ahead of the flame 


Si = Sf — u u (4.8) 

This gives a laminar flame speed of 1044 - 913 = 131 cm/s. 

An alternate method of finding the laminar flame speed is to apply conservation of mass to a 
control volume containing the flame. Transforming the moving flame to stationary coordinates by 
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Figure 4.3: Density profile, chemical energy release, and pressure across a one-dimension flame 
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subtracting the flame velocity, the continuity equation becomes 


Pb(Sp ~ life) - Pu(Sf ~ u u ) 


(4-9) 


Using equation (4.8) to replace Sp as *f u u , solving for the laminar flame speed, Si, gives 


Si — 


lit, - u b 


(4.10) 


(Pu/Pb) - 1 

Since the pressure is nearly constant, and the temperature ratio across the flame is 8.01, the density 
ratio pu/Pb , is also 8.01 from the equation of state (2.14). Then from equation (4.10), the laminar 
flame speed is calculated as 130.2 cm/s, a value very close to the computational value of 131 cm/s. 

The laminar flame theory, described in the previous section, provides an equation to calculate 
the laminar flame speed. This equation is reproduced below. 


s ‘ = /SkS)®' (4u) 

The I term is found from equations (C.4) and (4.2). After substituting the density from the equation 
of state, (2.14), and the mass fraction from equation (4.4), the integral becomes 

1 = (i) ! V' 4 ' £ (4.12) 

The theory assumes a constant value of the thermal conductivity. The present computation 
allows the thermal conductivity to be a function of temperature, as given by equation (2.27). To 
evaluate equation (4.11) for the laminar flame speed, an average value of temperature was used to 
determine the thermal conductivity. Using an unburned gas temperature of 293 K, the theoretical 
value of the burned gas temperature is 2344 K, as determined from equation (4.7). This results in 
an average temperature of 1319 K, which is used to calculate the thermal conductivity. Evaluating 
the integral, I, numerically using a pressure of 100 torr, the theory gives a laminar flame speed of 
103 cm/s. This value is 27% off from the value determined computationally. Kuo 2 states that the 
theory does not give very accurate results, but does predict the general trends that flame speed is 
proportional to the square root of the thermal conductivity. 

The laminar flame thickness is another important property of the flame, defined in equation (4.5) 
and shown again as 


x T b — T u 

* ~ (dT/dx) max (413) 

Evaluating this expression from the temperature plot gives a flame thickness of 0.023 cm. For the 
cell size of 0.001 cm, this gives 23 cells in the flame thickness. From laminar flame theory, equation 
(4.6) provides a formula for calculating the thickness. This formula relates the flame thickness to 
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Table 4.1: One-dimensional Flame Properties 



Si (cm/sec) 

6 t (cm) 

T b /T 0 

CMRFAST2D 

131 

0.023 

8.01 

Khokhlov et al. 37 

144 

0.022 

8.00 

Theory 67 

103 

0.024 

8.00 

Experimental 37, 66 

150 

~0.020 

NA 


the thermal conductivity and the flame speed. Using the values from the calculation of the flame 
speed, the theoretical value of the flame thickness is 0.024 cm, a 4.2% difference from the computed 
value. 

Table 4.1 compares the values of the flame speed, thickness, and temperature ratio obtained from 
the computed results using CMRFAST2D, separately computed results of Khokhlov et al., 37 the 
experimental values 37,66 which the chemistry model was based on, and the theoretical results. The 
two computed results agree well with each other and with the experimental results. The theoretical 
results are the furthest from the three other results, but as stated above from Kuo, 2 the theory 
provides the trends, but is not very accurate. 

4.4 Transport Property Test 

These tests compare computations of the laminar flame with the transport properties. p r , k r , 
and D ry halved and doubled, so that the Lewis and Prandtl numbers were always unity. The goal of 
these tests is to show the response of the flame to a change in the thermal conductivity. Figure 4.6 
shows the temperature, velocity, and mass fraction profiles for the two cases. Shown for comparison 
is the results for the standard case flame. The x axis was repositioned so that all three cases could 
be plotted in the same figure. This was done by defining x = 0 as the point where the maximum 
value of the temperature derivative, ( dT/dx) max occurred. 

The temperature profiles shows that the ratio of burned to unburned gas temperature is essen- 
tially unaffected by the thermal conductivity. This is because the ratio is a function of chemical 
parameters, as shown in equation (4.3). The ratio is independent of the thermal conductivity. This 
is also true of the mass fraction profiles. 

The velocity profiles show that when the thermal conductivity is doubled the velocity increases 
and when the thermal conductivity is halved the velocity decreases as compared to the standard 
case described in the previous section. To obtain the flame speed, the continuity method as given by 
equation (4.10) can be used. These values along with the flame thickness as determined by equation 


46 





Figure 4.6: Temperature, velocity, and fuel mass fraction for the cases of the thermal conductivity 
halved, k/2, and doubled, 2k. 
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Table 4.2: Laminar Flame Speed and Thickness Variations with Thermal Conductivity 



k 

2k 

\/2(k value) 

k/2 

(k value) / y/2 

St (cm/s) 

131 

194.4 

185.3 

90.7 

92.6 

S t (cm) 

0.023 

0.032 

0.033 

0.016 

0.016 


Table 4.3: Grid Study Flame Thickness and Velocity 


cell size (cm) 

St (cm/s) 

S t (cm) 

0.0010 

131 

0.023 

0.0015 

131 

0.023 

0.0020 

128 

0.025 

0.0030 

125 

0.030 

0.0040 

125 

0.036 


(4.5) are shown in Table 4.2. 

The laminar flame theory described in the previous section shows that the flame thickness and 
speed are directly proportional to the square root of the thermal conductivity. From equations (4.1) 
and (4.6), this is shown as 

Si oc Vk (4-14) 

S t oc y/k (4.15) 

If the thermal conductivity is doubled, 2k, the laminar flame speed should increase by (%/2)(131) = 
185.3, which is within 4.9% of the calculated value. Similarly, when the thermal conductivity is 
halved, k/2, the flame speed should decrease to (131)/(\/2) = 92.6, which is within 2.1% of the 
calculated value. The same trends are true for the flame thickness which for both cases are nearly 
identical to the calculated values. 

4.5 Grid Study 

For these test cases the cell size was increased to study the effect of the gridding on the results 
obtained above. Figure 4.7 shows the temperature, velocity, and mass fraction for five cell sizes 
ranging from 0.001 cm to 0.004 cm. Table 4.3 gives the calculated flame thicknesses and the laminar 
flame speed. The maximum error in flame speed is 4.6%. The errors for the flame thickness range 
from 8.7% at 0.002 cm to 56% at 0.004 cm. 

The simulations for the two-dimensional channel flame use a cell size of 0.002 cm in the lengthwise 
direction. This allowed about 13 cells in the flame thickness area. The decreased accuracy for this 


48 



cell size, 2.3% in flame speed and 8.7% in thickness, is offset by the increased step sized and domain 
size increases for the same number of cells. 
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Figure 4.7: Temperature, velocity, and fuel mass fraction varying grid resolutions 
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Chapter 5 


Two-Dimensional Channel Results 


This chapter describes the results for the two-dimension channel configuration. The wall bound- 
ary conditions for both an adiabatic and an isothermal wall will be discussed. For all cases, the 
boundary condition for the species at the wall is noncatalytic. 


5.1 Adiabatic Wall : Planar Ignition 


For these studies, the boundary condition for the temperature at the waill is adiabatic. The 
ignition method is described in section 2.8 : a discontinuity in temperature, species mass fraction, 
and density located 0.3 cm from the end wall of the tube. Table 5.1 lists the flame simulations that 
were performed. In the table, h is the channel height, /i r and D r are the reference values for the 
transport properties from equations 2.26 and 2.29, Re is the Reynolds number, and Le is the Lewis 
number. The reference Reynolds number is based on the laminar flatme speed, channel height, and 
kinematic viscosity in the unburned gas. For all the simulations the reference value for the thermal 
conductivity, fc r , is kept constant. 


Table 5.1: Adiabatic Wall Studies : Planar Ignition 


Case 

h (cm) 

fi r X 10 -8 

D T x 10- 8 

Re 

Le 


0.254 

2.439 

2.439 

76.2 

1.0 

AD2 

0.254 

4.878 

2.439 

38.1 

1.0 

AD3 

0.508 

2.439 

2.439 

152.4 

1.0 

AD4 

0.254 

2.439 

1.951 

85.2 

0.8 

AD5 

0.254 

2.439 

2.927 

69.6 

1.2 
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5.1.1 Case ADI : Baseline 


This first case is the baseline case that will be used for comparisons. The channel aspect ratio, 
length to height, is 32.2. Figure 5.1 shows the temperature profiles for the first 0.6 ms. The initial 
condition is a discontinuity of temperature, density, and mass fraction, as described in section 2.8. 
From this discontinuity in temperature, a flame develops. As the flame forms and begins to propagate 
down the channel, the wall inhibits the flame motion. The result is that the flame at the centerline 
of the channel moves faster than near the wall, resulting in the structure shown in Figure 5.1. In 
time, the curved portion of the flame near the wall grows longer. By 0.607 ms, this curved portion 
is about 0.3 cm long. 

Figure 5.2 shows the flame position as a function of time at two positions in the tube: the 
centerline and the wall. The position of the one-dimensional flame (shown in Figure 4.4) is also 
shown for comparison. The flame position is defined as the location in the flame front where the 
temperature is 1000 K. The flame position at the centerline closely follows the position for the one- 
dimensional flame for approximately 0.3 ms and then it becomes larger. The position of the flame 
at the wall is initially less then the one-dimensional result, but by about 0.55 ms it also is greater 
than the one-dimensional location. Figure 5.3 shows the flame velocity relative to the laboratory 
for the flame at the centerline and the wall, and these are compared to the velocity for the one- 
dimensional flame (also given in Figure 4.5). The flame speed has greatly increased over that of 
the one-dimensional flame. The one-dimensional flame reaches a constant velocity, but both the 
two-dimensional centerline and wall flames are still increasing at the end of the channel. 

Figures 5.4 - 5.11 show the temperature (T), density (p), vorticity (u>), pressure (p), velocity in 
the x direction (u), velocity in the y direction (v), velocity vectors, and instantaneous streamline 
pattern superimposed on the temperature contours at times of 0.433 ms, 0.758 ms, 1.104 ms, and 
1.450 ms. The vorticity is defined as twice the angular velocity and is given by 


dv du 
dx dy 


(5.1) 


These results are shown on the same scale to facilitate comparisons at the different times. 

The temperature profiles maintain the curved shape shown in Figure 5.1. By 1.450 ms, the 
curved flame has grown to be nearly 1.5 cm long. The flame thickness in the curved region and near 
the centerline are approximately the same, they only appear different due to the different scaling of 
the x and y axes. The shape of the density contours closely follows the shape of the temperature 
contours, because the pressure is nearly constant through the flame. The mass fraction profiles are 
not shown since, for adiabatic conditions, those profiles also resemble the the temperature profiles. 
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This can be shown by solving equation (4.4) for the mass fraction Y y and substituting in equation 
(4.3) to find 


Y = 


T h -T 


(5.2) 


T b -T u 

On Figure 5.4, the pressure slightly increases near 2.5 cm and this extends to about 3 cm. This 
is due to a reflected wave caused by the initial pressure wave reaching the outflow boundary and 
reflecting back into the domain. 

The u velocity shows that the speed of the gas in the unburned region increases in time. The 
v velocity shows that there is an upwelling at the wall in the curved-flame region. The maximum 
v velocity in this upwelling region remains essentially constant as the flame moves down the entire 
channel. The explanation for this upwelling lies in the observation that the material expands as it 
burns. The burned material cannot move downward because of the solid wall, and it cannot move 
backward, as it is being blocked by the previous burned material. The only direction for the material 
to expand is forward toward the outflow end of the channel. The burned material near the lower wall 
is restricted in moving forward because of the boundary layer ahead of the flame. The material that 
burned near the lower wall that is unable to move forward must then expand upwards toward the 
centerline. The material flows upwards behind the the flame, turns at the centerline, and provides 
additional momentum to the center of the flame. 

The instantaneous streamlines, or fluid paths, are plotted with the temperature contours to 
show the streamlines relative to the flame location. These figures show the flow upwelling and then 
turning toward the outflow boundary. The highest values of vorticity occur in the boundary layer. 
There is relatively little vorticity behind the flame and near the centerline of the channel, most of 
the vorticity occurs in the boundary layer and along the curved flame surface. 

Figure 5.12 shows the evolution of the boundary-layer thickness along the channel in front of 
the flame. The times in the figure are about 0.2 ms apart. The boundary layer thickness, <5 99 , is 
defined as the vertical height above the wall where the u velocity reaches 99% of the centerline value. 
At a given location, the boundary- layer grows in time until the flame reaches that location. As the 
flame moves over that location it destroys the boundary layer. Initially the boundary-layer thickness 
grows rapidly and then this growth rate slows. The centerline velocity, Figure 5.14, shows that the 
fluid velocity in the unburned gas accelerates as the flame moves down the channel. At 0.217 ms 
the centerline velocity drops to zero at about 7 cm. This is because this is the furthest point that 
the initial pressure wave, traveling at the local sound speed, has reached. Ahead of this location, at 
this time, the gas is unaware of the developing flame. 

Figure 5.13 shows the displacement thickness in front of the flame for various times. The dis- 
placement thickness represents the distance that the inviscid freestream is displaced from the wall 
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due to the decrease in velocity in the boundary layer. It is defined as 



This plot is similar to the boundary layer thickness but it shows that the displacement thickness 
levels off more quickly, to a value close to 0.017 cm. The small displacements in the profiles at 0.433 
ms and 0.650 ms represent the reflection of the initial pressure wave from the outflow boundary and 
the reflection of the reflected wave from the flame front. This is also seen in the centerline velocity 
discussed in the previous paragraph. 

Figures 5.15 and 5.16 show the boundary layer growth at different channel locations. The time, 
t x , is the local time at a channel location, x . This is defined as the time starting from when the 
initial pressure wave, traveling at the local speed of sound in the unburned gas, has reached that 
point. The local time, t x , is calculated from the globed time ,£, as 


(x - 0.3) 

\J~lR t T u 


(5.4) 


The 0.3 is needed since the initial high temperature region used to initialize the flame was located 
up to 0.3 cm, so this is where the initial pressure wave is formed. 

Examining these figures, the boundary layer growth curves at different channel locations collapse 
onto the same curve for times greater than about 0.3 ms, becoming self-similar. The solid line in 
both figures gives the theoretical results from Stokes’ First Problem discussed in section 3.2.3. From 
White, 53 the boundary layer thickness for this problem is given by 


$99 = 3.64 yjv u t x 


(5.5) 


In the theory, it was assumed that the freestream velocity and kinematic viscosity, v u , are constant. 
Using equation 2.26, the kinematic viscosity is given as 



Mr 

P 



(5.6) 


The theoretical value of the boundary layer thickness was calculated with the kinematic viscosity 
held constant at a value based on the unbumed gas conditions. 

The boundary layer growth initially follows the theory, but it then slows its growth. This is due 
to two effects acting together. The primary cause for the deviation is that the centerline velocity 
is not constant (Figure 5.14), but increases in time for all channel locations ahead of the flame. 
Figure 5.17 gives the centerline velocity ahead of the flame for the local time at selected channel 
locations. After about 0.2 ms, the curves for the centerline velocity have a similar shape, which gives 
similar accelerations at a given channel location at the same local time. The acceleration of the flow 
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causes the boundary layer thickness to decrease from that predicted by Stokes’ First Problem. This 
is shown in Appendix A for three different accelerations. It is the acceleration of the centerline 
velocity that is the primary cause of the slower boundary layer growth. 

The secondary cause of the slower growth is due to the change in the kinematic viscosity in the 
unburned gas. For the theory, it was assumed that the kinematic viscosity was constant at the initial 
conditions of the unbumed gas. The pressure waves that are generated as the flame propagates down 
the channel cause a decrease in the value of the kinematic viscosity in the unburned gas. The pressure 
wave increases the temperature and density above the initial values. The relative increase in density 
is greater than the relative increase in temperature (shown in Appendix B). By equation 5.6, the 
greater rise in density lowers the kinematic viscosity. Since the boundary layer growth is proportional 
the square root of the kinematic viscosity, a lower value gives a slower boundary layer growth. 

5.1*2 Case AD 2 : Effects of Increased Viscosity 

This second case has initial conditions identical to case ADI described in the previous subsection. 
The aspect ratio of the channel is the same, 32.2. The reference coefficient for the viscosity is 
doubled, since the viscosity controls the growth rate of the boundary layer. By doubling the reference 
coefficient, the Reynolds number based on laminar flame speed is halved. 

The development of the flame from the initial conditions is similar to the development in the first 
case, Figure 5.1. The flame forms from the discontinuity in temperature and starts to propagate 
down the channel. The wall inhibits this motion and a curved flame develops. 

Figure 5.18 shows the flame position as a function of time for the channel centerline and wall. 
Figure 5.19 shows the velocity of the flame relative to laboratory coordinates. From these two figures 
it can be seen that the flame is faster then in the case ADI and therefore exits the channel earlier. 
By 1 ms, the flame has reached 1 cm further down the channel and is moving 3000 cm/s faster than 
the first case. 

Figures 5.21 and 5.22 show the temperature, density, pressure, vorticity, velocities, and instan- 
taneous streamlines at a time of 1.100 ms. This time is very close to the time of 1.104 ms in 
figures 5.8 and 5.9 for case ADI. Both sets of figures are also in the same scale to facilitate com- 
parisons. Examining these two sets of figures shows that the flame has moved further down the 
channel then in the first case. Comparing the temperature plots, the length of the curved region 
of the flame near the wall has increased for the second case. Figure 5.20 shows the length of this 
curved region for the two cases, Axj. This length increases faster for the doubled viscosity case. 
The added length provides more burned material in the fountain effect, which in turn gives a greater 
push to the flame. This is what accelerates the flame faster then in the first case which has less of 
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a curved flame region at the wall. 

Figure 5.23 shows the evolution of the boundary layer thickness in front of the flame. The times 
shown are similar to the times in Figure 5.12 for comparisons. Comparing these two figures it can 
be seen that the boundary growth is greater for the higher viscosity case, as is to be expected. 
Figure 5.24 shows the time history of the displacement thickness along the channel. The plot 
is qualitatively similar to Figure 5.13 in case ADI, but the magnitudes are greater. Both the 
displacement and boundary layer thicknesses have grown thicker in less time, providing a greater 
blockage to the flame induced flow, which also adds to the increased fluid velocity ahead of the flame. 

Figure 5.25 and Figure 5.26 show the boundary layer growth at different channel locations at 
the local time, t x . The boundary layer growth is again seen to become self-similar in local time but 
the deviation from theory is greater then was shown for case ADI. This greater deviation is caused 
by the greater acceleration of the flow. Figure 5.27 shows the velocity in the channel for different 
times. Comparing this to case ADI, Figure 5.14, the velocity is seen to be faster at a given time. 
Figure 5.28 shows the velocity at different channel locations for local time. The velocity curves are 
similar in shape after about 0.2 ms, which leads to similar accelerations at a given channel location 
for the same local time. 

5.1.3 Case AD3 : Effects of Increased Channel Height 

For this case, the channel height is doubled, to 0.508 cm. This has the effect of doubling the 
Reynolds number, to 152.4, and halving the channel aspect ratio, to 16.1. The initial conditions 
are identical to case ADI, described in section 5.1.1. The reference coefficients for the transport 
properties are the same as case ADI. 

Figure 5.29 shows the temperature profiles for the first 0.6 ms from a flame forming at the initial 
temperature discontinuity located at 0.3 cm. The times shown are similar to the times in figure 5.1 
for the first case. Comparing the two figure shows that the flame in AD3 has not moved as far down 
the channel in the same amount of time. Figure 5.30 shows the position of the flame versus time. 
The flame position is defined as the location of the 1000 K isotherm. Now the flame takes about 2.6 
ms to move out the end of the channel, as opposed to about 1.5 ms for the first case. Figure 5.31 
shows the velocity of the flame relative to laboratory coordinates. The flame is slower moving than 
case ADI, which accounts for the longer time to reach the end of the channel. 

A second interesting feature shown in figure 5.29 is the formation of a bulge at the flame front 
that begins near the lower wall and moves upward towards the channel centerline in time. The 
vertical location of the furthest downstream point of the flame front is shown in figure 5.32. By the 
end of the channel, the vertical height of the bulge has moved to about 0.165 cm, or a little over half 
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the distance to the centerline. This bulge creates a flame that bends back at the centerline, forming 
a cusp. Since the bulge is moving upwards, the cusp forms, grows, and then begins to disappear. 
This is shown in figure 5.33, which gives the horizontal distances between the position of the flame 
front at the centerline to the position of the flame front at the wall, the position at the the bulge 
to the wall, and the position at the bulge to the centerline. For comparison purposes, the distance 
from the centerline to the wall for the first case, ADI, is also shown. The cusp, or the bulge to 
centerline distance, forms, grows to a maximum at about 1 ms, and then begins to disappear. This 
figure also shows that the curved flame front near the wall (bulge to wall distance) is not as long as 
in case ADI. Since this region is where the upwelling originates from, a smaller region provides less 
upwelling which provides less additional momentum to the flame. 

Figures 5.34 - 5.41 show the flowfield properties at times of 0.541 ms, 1.061 ms, 1.601 ms, 
and 2.121 ms. All the figures are plotted in the same scale for comparison purposes. Due to the 
increase domain height though, the scale is not the same as the related figure from case ADI. The 
temperature and the density plots show the formation of the bulge and of the cusp at the centerline. 
The bulge moves upwards in time as the cusp grows and then shrinks in size. The pressure plots 
show the initial waves moving out from the flame. The pressure change at later times is about 0.1%. 
The vorticity shows that the highest values occur in the boundary layer and along the flame surface, 
which is where the flow curvature is the highest. 

The velocity plots show how the bulge forms. In figure 5.35 the curved flame front near the wall 
is short. There is therefore a small amount of upwelling. This small amount of upwelling provides 
additional momentum to the flame that is closer to the wall then near the centerline. The flame at 
the centerline receives less additional momentum and moves more slowly then the flame closer to 
the wall, producing a bulge. As the flame propagates, the curved flame region near the wall grows 
longer and it produces more upwelling which is distributed over more of the flame surface. The 
bulge moves upwards as the push is felt over a greater area of the flame surface. The cusp at the 
centerline begins to shrink. As time increases, the upwelling reaches the centerline and is distributed 
over the entire flame surface. By the end of the channel, the bulge is about 0.035 cm ahead of the 
centerline, about 1^ laminar flame thicknesses. 

Figure 5.42 shows the time evolution of the boundary layer thickness along the channel ahead 
of the flame. Figure 5.43 shows the displacement thickness along the channel. Compared to the 
thicknesses in case ADI (figures 5.12, 5.13), the boundary layer has grown thicker for a given time. 
Figure 5.44 shows the centerline velocity profiles. Compared to case ADI, (figure 5.14), the velocity 
is less. The smaller foot causes less of an upwelling which gets distributed over the larger flame area. 
This results in less of a momentum push to the flame which gives a lower acceleration and therefore 
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a slower velocity. With the lower acceleration, the boundary layer can grow thicker compared to the 
first case. A secondary effect which acts to lower the acceleration comes from the blockage effect 
of the boundary layer. As the boundary layer grows, the effective channel height decreases. Even 
though the boundary layer is thicker, the relative decrease in the channel height is less due the 
increased channel height. 

Figures 5.45 and 5.46 shows the boundary layer growth related to the local time. The local time 
was defined in equation 5.4 as the time starting from when the initial pressure wave first reaches 
that point of the channel. As with the previous two cases, the boundary layer growth becomes self- 
similar in local time. The comparison to the solution of Stokes’ First Problem, <$99 = 3.64 yjv u t x , is 
better then the previous cases. This is due again to the lower acceleration. Stokes’ solution assumed 
zero acceleration. This case has an acceleration that was lower then the previous two cases, so the 
comparison is better. 

5.1.4 Case AD4 : Effects of Decreased Mass Diffusion 

This case has a channel geometry the same as case ADI. The channel aspect ratio is 32.2. This 
case has the reference coefficient for the mass diffusion decreased so that the Lewis number decreases 
to 0.8. The Lewis number was defined as the ratio of the mass diffusivity to the thermal diffusivity, 
and was given in equation 2.30. The Lewis number was decreased, since according to Glass man , 5 7 
flames instabilities may occur for Lewis number not equal to unity. From laminar flame theory, the 
laminar flame speed is related to the Lewis number by 



Compared to the laminar flame speed from the first case, ADI, this flame will have a higher flame 
speed, about 1.12(Si)adi- For this case no flame instabilities were present even though the Lewis 
number was less then one. 

The initial development of this flame from the planar discontinuity is very similar to the first 
case, as shown in figure 5.1. As with the first case, no cusp forms at the centerline as occurred in case 
AD 3 . Figures 5.47 and 5.48 shows the flame position and velocity relative to laboratory coordinates. 
The flame exits the channel at a slightly earlier time as compared to the first case, about 0.02 ms 
earlier, and has accelerated to a slightly higher speed, about 160 cm/s higher. 

Figures 5.49 and 5.50 show the temperature, density, pressure, vorticity, velocity components, 
and instantaneous streamlines at a time of 1.172 ms. The figures show that the flame resembles the 
basic structure that was seen in all the previous cases. The flame is curved near the lower wall which 
produces the upwelling, which drives the flame faster down the channel. The pressure waves that 
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are generated by the flame are planar and move out ahead of the flame toward the outflow. The 
vorticity is concentrated in the boundary layer and along the flame surface, since these are where 
there is flow curvature. 

Figure 5.51 shows the evolution of the displacement thickness along the channel. The times are 
similar to the times in figure 5.13 for comparisons. Figure 5.52 shows the boundary layer thickness 
along the channel and figure 5.53 shows the how the centerline velocity increases in time as the flame 
accelerates down the channel. These three figures are all qualitatively similar to the boundary layer 
profile plots of the first case. Figures 5.54 and 5.55 give the boundary layer growth at different x 
locations for varying local time. As with the other cases, the results are self-similar due to the similar 
slopes of the centerline velocities in local time as shown in figure 5.56. Due to the acceleration, the 
boundary layer growth is less than predicted by Stoke’s First Problem. 

5*1.5 Case AD5 : Effects of Increased Mass Diffusion 

This case again has a channel geometry the same as case ADI, with a channel aspect ratio 
of 32.2. The reference coefficient for the mass diffusion was increased so that the Lewis number 
increased, Le > 1. Compared to the laminar flame speed from the first case, this flame will have a 
lower speed, about 0.91(S/)^pi. The initial development of the flame from the planar discontinuity 
is very similar to the first case, ADI. As with the case in the previous section, no centerline cusp 
appears in the flame. Figure 5.57 shows the flame position in time. Figure 5.58 gives the flame 
speed relative to laboratory coordinates. The flame exits the channel at a slightly lesser time then 
the first case and is traveling at a slightly lower speed. 

Figure 5.59 and 5.60 show r the flowfield properties for this flame at 1.284 ms. The basic structure 
of the flame resembles the flame in the previous section and for case AD 1 . The upwelling at the wall 
provides the added momentum push to the flame at the centerline which accelerates it down the 
channel. Figure 5.61 and 5.62 show the displacement thickness and boundary layer thickness growth 
in time along the channel. Again the results are similar in shape to the first case. Figure 5.63 and 
5.64 show the growth of the boundary in local time at selected channel locations. The curves for the 
different channel locations are self similar in local time and fall below Stokes’ First Problem due to 
the flow acceleration. 

5.1.6 Comparisons of the Planar Ignition Cases 

Figure 5.65 shows the evolution of the relative flame area in time for the five planar ignition cases. 
The relative flame area is the flame surface area, defined as the length of the 1000 K isotherm, divided 
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by the channel cross-sectional area. Since the simulations were two-dimensional, a unit length into 
the page was assumed. Figure 5.66 shows the evolution of the flame speed relative to laboratory 
coordinates. Figure 5.67 gives the evolution of the boundary layer thickness at the 8 cm location. 
Also shown on this figure is the theoretical solution of Stokes’ First Problem. Comparisons will 
made relative to the first case, ADI. 

The second case, AD2, has the viscosity doubled. The increased viscosity causes the boundary 
layer to grow thicker for a given amount of time. This thicker boundary layer causes the curved 
area near the wall to be longer which leads to a larger flame area. The larger flame area burns 
more material which provides greater additional momentum. The greater additional momentum 
accelerates the flame faster which causes the boundary layer to grow slower then predicted by 
Stokes’ First Problem. 

The third case, AD3, has the channel height doubled. The relative burn area is lower at a given 
time then for any of the other cases. This lower burn area provides less of a push from the expanding 
burned material. This lower momentum push causes the flame to move slower. Since the flame has 
less acceleration, the boundary layer growth more closely matches that from Stokes’ Problem. Even 
though the boundary layer growth is greater in time then the first case, it has less of an effect on 
the flow due to the increased channel height. 

The fourth and fifth planar ignition cases, AD4 and ADS, have the channel height and viscosity 
the same as case ADI but the Lewis number has been changed. The fourth case has the Lewis 
number lowered and the fifth case has the Lewis number raised. For both of these cases, the burn 
area and flame speed are close to the first case, with AD4 being slightly higher and ADS being, for 
the most part, slightly lower. The higher flame speed for case AD4 leads to a higher acceleration 
which leads to it being slightly below the curve for ADI in figure 5.67. Case AD5, which has a 
lower acceleration is slightly above case ADI in that figure. Due to the acceleration though, both 
are below the curve for Stokes’ Problem. 

Figure 5.68 shows the flame speed plotted against the relative flame area. The relative burn area 
is defined as the ratio of the surface area of the flame front, Aj , to the channel cross-sectional area, 
A c . The curves for all cases fall nearly on top of each other. Using a control volume fixed to the 
moving flame, the continuity equation reduces to 

Pb(SF - Ub)A c = Pu(Sf - u u )Af (5.8) 

Assuming that the burned gas velocity is zero, using equation (4.8) for the definition of the laminar 
flame speed and the equation of state (2.14), equation (5.8) becomes 67 
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Table 5.2: Adiabatic Wall Studies : Spark Ignition 


Case 

h (cm) 

Re 

Le 

ASl 

0.254 

76.2 

1.0 

AS2 

0.508 

152.4 

1.0 


This equation states that the flame speed, relative to the laboratory, is directly related to the ratio 
of the flame area to channel cross-sectional area. In the above equation, both the temperature ratio 
and laminar flame speed are functions of the Lewis number. The relationship between the laminar 
flame speed and Lewis number was given by equation (5.7). According to Zel’dolvich, 67 an increase 
in Lewis number causes an increase in the temperature ratio and a decrease causes a drop in the 
ratio. 

For cases ADI - AD3, the Lewis number was unity. The temperature ratio is found from equa- 
tion (4.7) to be 8 from theory. Using the laminar flame speed of 128 cm/s from table 4.3, the 
theoretical slope for the first three cases is 1024. Performing a least squares fit, 68 assuming a zero 
intercept, gives a slope of 1083, a 5.8% difference from the theoretical value. Part of the difference 
could come from the error associated with the calculation of the speed from the position data. Also 
the theory assumes an inviscid flow*, but for all these cases the viscosity was non-zero. 

For case AD4 the Lewis number was decreased to 0.8. This has the effect of increasing the 
laminar flame speed and decreasing the temperature ratio. Using an estimation of the laminar flame 
speed of 143 cm/s, 1.12(5j)adi, and a temperature ratio of 7.8, the slope is 1116. Curve fitting the 
data, a slope is 1101, a difference of 1.4%. For case AD5 the Lewis number was raised to 1.2 and the 
effect was reversed, the laminar flame speed was increased and the temperature ratio was lowered. 
Using a laminar flame speed of 117 cm/s and a temperature ratio of 8.4, the theoretical slope is 982. 
From the curve fit the slope is 1069, a difference of 8.1%. 

5.2 Adiabatic Wall : Spark Ignition 

For all of these studies, the w’all boundary condition is adiabatic. The difference from the 
previously described cases is the ignition method. For these cases, the ignition is done with a 
simulated spark. A quarter-circular, high temperature region is placed at the center of the channel 
along the end wall. The reference coefficients for the transport properties, fi r , k r , D r , are all given 
in Table 2.1, so that the Lewis and Prandtl number for the all cases is unity. Table 5.2 lists the 
flame channel simulations that w r ere performed. 

To ignite a flame from a spark, a minimum amount of energy, must be supplied from the 
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spark. Prom Turns , 64 this is the amount of energy required to raise the temperature of the spherical 
spark to the adiabatic flame temperature. 

Ei = jnr 3 m p bCv (T b -T u ) (5.10) 

This requires knowledge of the minimum spark radius, r m . Turns 64 outlines an approach that equates 
the rate of chemical heat release in the spark volume to the rate of heat loss due to conduction over 
the spark surface and estimates the minimum radius to be 

r m = 1.22 St (5.11) 

Using a different approach, Glassman 57 equates the chemical reaction time to the conduction time 
and estimates the minimum radius to be 


r m = 1.85* (5.12) 

Using the l amin ar flame thickness of 0.023 cm from chapter 4, the above two equation give a minimum 
radius of 0.028 cm and 0.043 cm, respectively. The radius used for the spark simulations was 0.05 
cm. Also the temperature of the initial burned material inside the spark was raised to 2930 K. This 
was done to ensure that there was enough energy to initiate flame propagation. 

5.2.1 Case AS1 : Baseline 

This case has a channel geometry the same as the first planar ignition case, ADI. The channel 
aspect ratio, length to height ratio, is 32.2. Figure 5.69 shows the temperature profiles for the 
first 0.39 ms of the simulation. From the initial temperature discontinuity, the flame develops and 
spreads outwards in a radial fashion. The flame traveling down the end wall toward the bottom 
of the channel is slowed down and the flame develops a non-circular shape. By the time the end 
wall flame has reached the bottom of the channel, a distance of about 0.077 cm from the initial 
temperature discontinuity at the centerline, the centerline flame has moved out a distance of about 
0.473 cm. In time, the flame assumes the shape similar to the first case, ADI. 

The final time in Figure 5.69, 0.391 ms, corresponds to the fourth time, £ 4 , in Figure 5.1. In the 
same amount of time, the present flame has traveled further down the channel and has developed 
a greater flame area then the first case. This is shown in Figures 5.70, 5.71, and 5.72 which show 
the flame position, velocity, and area ratio in time. As shown in Figure 5.72, the relative flame area 
to the channel cross-sectional area is greater for this case, AS1, as compared to the first case, ADI. 
The greater flame area results in more burned material which in turn provides greater additional 
momentum to accelerate the flame. Figure 5.72 shows a constant relative flame area from about 
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0.3 ms to 0.5 ms. During this time, the flame receives a constant push and moves at a constant 
speed. When the flame area begins to grow again, the flame accelerates due to the increase in burned 
material. 

Figures 5.73 - 5.78 show the flowfield properties at 0.102 ms, 0.209 ms, and 0.305 ms. Fig- 
ures 5.73 and 5.74 show the flowfield for the flame still traveling down the end wall of the channel. 
The velocity field induced by the flame shows that the expanding burned gas flows both behind the 
flame, up to the centerline, and then out, as well as down and along the bottom wall. The streamline 
patterns in the last plot show this. The vorticity plot shows a boundary layer that starts to grow 
from the end comer of the channel due to the initial expanding gas flowing down along the end wall. 
Figures 5.75 and 5.76 show the flowfield when the end wall flame has reached the bottom of the 
channel. The expanding burned material is then forced to flow up and out along the centerline due 
to the bottom wall boundary layer, as in the planar cases. By 0.305 ms, the flowfield has become 
qualitatively similar to the flowfield that developed in the planar ignition cases. 

The pressure plots show a difference from the planar ignition cases. Since the flame is initialized 
as a circular discontinuity, the initial pressure wave is circular and not planar. This circular waves 
moves outward and is repeatedly reflected off of the channel side walls. The reflected waves also 
interact with each other inside of the channel. Figure 5.79 shows the pressure ahead of the flame 
normalized by the initial pressure at a time of 0.166 ms, a time about half way between the first 
two flowfield times. This plot shows that at a given channel location, the pressure varies across the 
channel. This is due to the the pattern of the reflected waves in the channel. Figure 5.74 shows that 
the streamline pattern oscillates due to the pressure oscillations. As the flame approaches a profile 
similar to case ADI, the pressure waves straighten out and this oscillating pattern flows out of the 
system. 

Figure 5.80 and 5.81 show the evolution of the boundary layer thicknesses. At earlier times, 
both boundary layer thicknesses oscillate, which is caused by the fluctuating velocities due to the 
pressure wave interactions in the channel. As time increases, these waves leave the system and the 
thicknesses smooth. Figure 5.82 shows the velocity along the channel centerline for various times. 
At earlier times, the velocity is constant for a distance ahead of the flame, followed by a decrease in 
velocity to the end of the channel. The first time plotted corresponds to when the relative flame area 
was nearly constant. As time increases, the area ratio increases, and the flame velocity increases. 

Figures 5.83 and 5.84 show r the boundary layer growth in local time. Since the initial pressure 
wave is curved, the starting point for the local time calculation is when the pressure wave along 
the wall first reaches a given location. As can be seen in both figures, the boundary layer growth 
is less then Stoke’s First Problem, due to the flow acceleration. Also seen is the fluctuations in the 
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boundary layer thickness at early times due to the pressure wave reflections. By about 0.25 ms local 
time, these waves have passed by all locations and the thickness smooth and become self-similar. 

5.2,2 Case AS2 : Effects of Increased Channel Height 

This case has the same channel geometry as the third planar ignition case, AD3, in that the 
channel height is doubled, making the channel aspect ratio 16.1. Figure 5.85 shows the initial flame 
development for the first 0.672 ms. The circular spark deforms greatly in time. The end wall of 
the channel acts like the side wall for the portion of the flame that is moving along it. As seen at 
time t 3 , the portion of the flame moving along the end wall has obtained a shape similar to the 
flame that moves along the channel side wall for later times, for example ts- As the end wall flame 
approaches the side wall, the portion of the flame traveling along the end wall overtakes the flame 
a short distance from the end wall. This creates a flame remnant, as seen at time £ 5 . This remnant 
becomes detached from the main body of the flame which propagates down the channel. In time, 
the remnant burns itself out. 

The reason that the remnant was not seen for the first spark case is that there was not enough 
time, or distance traveled along the end wall, for it to develop. As seen in figure 5.69 for the first 
spark case, the flame moving down the end wall reaches the side wall by about 0.230 ms. For case 
AS 1 , the flame moving along the end wall is nearly straight, as opposed to the highly curved flame 
that exists in the present case. When the straight flame nears the end wall, the flame along the wall 
accelerates to the corner as seen at time t 3 in figure 5.69, so that no remnant forms. 

Figures 5.86 and 5.87 show the flame position and velocity. Compared to the third planar ignition 
case, AD3, which has the same channel height, this case exits the channel sooner and reaches a higher 
flame speed. This is caused by a higher burn area, as shown in figure 5.88. This figure compares 
the relative bum area for the two cases. The higher bum area for the spark case leads to a higher 
flame speed. Also, comparing figures 5.87 and 5.88, the flame speed has the same general shape as 
the area ratio. 

Figures 5.89 - 5.96 show the flowfield properties at 0.348 ms, 0.500 ms, 0.650 ms, and 0.801 
ms. Figure 5.89 and 5.90 show the flame when it is still moving down the end wall. The end wall 
flame has the same shape as the side wall flame seen in the previous cases. The velocity vectors 
and instantaneous streamlines show that the burned material from the end wall flame expands both 
ahead of and behind the flame. The flow that expands ahead of the flame turns the comer and 
proceeds out the channel. The material that expands upwards along the end wall then turns and 
provides an extra push to the flame near the centerline. This flow pattern is similar to the flowfield 
reported by Gonzalez 16 for the flow induced during a spark ignition. The vorticity plot shows that 
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the boundary layer for this case begins at the corner of the channel. The pressure plot shows that, 
like the first spark case ASl, the pressure contours are not straight. This is due to the initially 
circular pressure wave reflecting off of the side walls. 

Figures 5.91 and 5.92 shows the end wall flame as it nears the side wall. The expanding gas that 
was pushed ahead of the flame and around the corner is now pushed either along the flame surface 
or is expanded behind the flame. With the gas flow ahead of the end wall flame stopped, the flame 
traveling adjacent to the wall moves into a nearly stagnant gas. The boundary layer along the end 
wall formed by the expanding gas disappears when the velocity goes to zero. 

Figures 5.93 and 5.94 show how the remnant flame forms. The flame has reached the side wall 
of the channel. The velocity vectors show that the remnant is drawing material upwards towards it. 
This upward flow provides an added momentum push from the remnant flame towards the centerline 
of the channel. Due to the presence of the end wall, the additional push on the flame attached to 
the end wall is less due to the no-slip boundary condition bringing the velocity to zero. This allows 
the attached flame to overtake the flame in the center of the remnant. This is similar to what 
happens to a flame in a closed channel when it reaches the wall opposite from the ignition source. 
The flame that reaches the side wall is quenched, and a remnant is left behind by the main flame as 
it propagates down the channel. Figures 5.95 and 5.96 show the flame after the remnant has burned 
out. The flame assumes a shape similar to the shape in case ADI. Unlike case AD3, no cusp forms 
at the centerline. This is because the centerline flame always has an additional momentum push 
and can remain ahead of the rest of the flame body. 

The evolution of the displacement thickness and boundary layer thickness is shown in fig- 
ures 5.98 and 5.97. The oscillations at the earliest time shown are due to the pressure oscillations 
caused by the reflecting waves from the initial circular pressure wave. This is the same thing that 
happened in the first spark case, AS1. By later times the fluctuations have exited the channel and 
the thicknesses smooth out. Both thicknesses are greater than the first spark case due to the lower 
velocity. They are thinner though than the third planar case, AD3, because compared to that case 
this case has greater acceleration. The greater acceleration thins out the boundary layer. 

Figure 5.99 and 5.100 give the growth of the boundary layer in local time. Due to the acceleration 
of the flow, the boundary layer growth is less than Stokes’ First Problem until about a local time 
of 1.2 ms when the growth becomes greater. This corresponds to the deceleration of the flame as 
shown in figure 5.87. The oscillations that occur in the beginning of the plots are from the pressure 
reflections as described in the preceding paragraph. 
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5.2.3 Comparison of the Spark Ignition Cases 


Figure 5.101 show the evolution of the area ratio for the two spark ignition cases. Figure 5.102 
shows the time history of the flame speed. On both figures, the related case from the planar ignition 
simulations with the same channel height is also shown. For the first spark case, AS1, the area 
starts out less then the area for the planar ignition, due to the smaller spark size, but very rapidly 
surpasses the area from the first planar ignition case. Since the burn area is increased, the speed of 
the flame is increased. The second spark case, AS2, also starts out with less flame area, but quickly 
surpasses the bum area from the third plan air case. It always remains above the area for the third 
planar case. Figure 5.103 shows the growth of the boundary layer at the 8 cm location for both 
spark cases. Both are seen to be less then the solution to Stokes’ First Problem, except for the later 
part of AS2 which starts to grow faster. This is due to the deceleration of the flame as shown in 
figure 5.102. 

Figure 5.104 shows the flame speed and bum area plotted on the same graph. The curve for the 
flame speed is seen to follow the curve for the bum area. Figure 5.105 shows the flame speed plotted 
against the bum area. In both cases, for the initial expansion of the spark, the curves are nearly 
linear and lie almost on top of each other. The relationship between the flame speed and bum area 
was given in equation 5.9, and is repeated here. 

(l)K£) <5i3) 

Using a laminar flame speed of 128 cm/s and a temperature ratio of 8, the theoretical value for the 
slope is 1024. For the initially linear portion of both curves, the slopes are 991 for ASl and 970 
for AS2, differences from theory of 3.2% and 5.3% respectively. The curve for ASl, after the initial 
expansion of the spark, the bum area curve flattens out, as does the flame speed. Since the two 
curves are slightly out of phase with each other this accounts for the notch back at an area ratio 
of about 5. After this constant region the flame area begins to increase again with a slope of 1096, 
7.0% from the theoretical value. After the spark expansion, the bum area decreases, then slightly 
increases again. Since there is a lag time in the response of the flame speed, this accounts for the 
wandering of the curve in figure 5.105 after the initial straight section. 

5.3 Isothermal Wall 

This case has a channel geometry the same as the first adiabatic wall case, ADI. The channel is 
0.254 cm in height and 8.2 cm in length. The ignition method is a planar high temperature region 
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located at the closed end of the channel from 0 cm to 0.3 cm. This is the same ignition method 
as in the first adiabatic wall case. The wall temperature boundary condition from 0.3 cm to the 
outflow end of the channel is isothermal, with a wall temperature of 293 K, which is equal to the 
initial unburned gas temperature. The species boundary condition is non-cat alytic. The reference 
coefficients for the transport properties, p r , kr, and D r , are given in Table 2.1, so that the Lewis 
and Prandtl numbers are unity. 

Figure 5.106 shows the flame position in time at the centerline of the channel. Shown for 
comparison is the flame position for the first adiabatic wadi case and the 1-D flame. Figure 5.107 
gives the isothermal flame speed in time. Compared to case ADI, the isothermal flame takes a longer 
time to exit the channel. The flame speed also shows that the flame accelerates and decelerates as 
it propagates down the channel. This is consistent with the “jerky” motion of a flame, 9,27 where 
the motion oscillates. It is also interesting to note that the mean propagation location and velocity 
is close to the 1-D flame results. 


Figure 5.108 shows the flowfield temperature, T, fuel mass fraction, F, chemical energy release, 
A£ c , and pressure, p, at 0.217 ms. The chemical energy release is the amount of energy that is 
transferred from the zero point energy term, pqY, to the sensible energy term, pCyT in the energy 
equation, (2.25), during the chemical reaction. 


ait dY 

A£ t = 


(5.14) 


The pressure waves that emanate from the initial planar discontinuity are nearly planar, as with the 
first adiabatic wall case. As seen in the temperature plot, the cold wall cools the hot burned product 
behind the flame. The mass fraction profiles do not follow the temperature profiles as they did for 
all of the adiabatic wall cases. This is due to the mass diffusion effects. The energy release shows 
that the flame is quenched, or extinguished, due to the cold wall. This is evident by the energy 
release not touching the lower wall. The quenching distance for a flame traveling between parallel 
plates is derived in Appendix C and is repeated here. 


d 9 = 6 t j2l^. 


'T q - T u 
T b -T u _ 


where d q is the quenching distance, and T q is a quenching temperature. Using a value of 1000 K 
for the quenching temperature, and evaluating the thermal conductivity, k , at the average flame 
temperature of 1319 K, the theoretical quenching distance is 0.011 cm. Again using the 1000 K 
isotherm as a reference, the average separation distance over the entire calculation between the 
wall and the nearest position of the isotherm to the wall is 0.005 cm. Since this value represents 
quenching for half of the flame, due to the symmetry condition at the centerline, the total quenching 
distance is 0.010 cm, about a 10% difference from the theoretical value. 
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Figure 5.109 shows the velocities and instantaneous streamlines for this time. The expanding 
burned gas induces a flow ahead of the flame in the unburned material toward the open end of the 
channel. A flow also develops behind the flame in the burned region. This is due to the cooling of 
the hot burned product by the cold wall. This cooling of the gas by the wall causes the density to 
decrease. This decrease in density provides space into which the burned material can expand. 

Figures 5.110 and 5.111 show the flame properties at 1.304 ms. By this time in the adiabatic wall 
cases, the flame is nearly at or has exited the channel. For this case, the flame has reached about 
1.5 cm down the channel from the initial planar discontinuity. The temperature plots shows that 
the cooling is more pronounced, with a large cool region existing behind the flame. This large cool 
region behind the flame allows the burned material to expand back toward the end wall. Examining 
the u velocity plot, there is a greater magnitude of reverse flow then there is for flow out of the 
channel ahead of the flame. The streamlines show a region of reversed flow has developed ahead of 
the flame close to the lower wall and extends to the exit of the channel. The flow in the center of 
the channel is still directed out of the channel. The pressure plots shows that there is an adverse 
pressure gradient in the channel. 

Figures 5.112 and 5.113 show the flame properties at 1.760 ms. The streamlines show that the 
reversed flow region near the lower wall is gone. The flow ahead of the flame is entirely moving 
toward the exit of the channel. The velocity at the centerline ahead of the flame has accelerated 
to a faster speed. The pressure gradient ahead of the flame has also reversed itself into a favorable 
one. Figures 5.114 and 5.115 give the flame properties at 2.519 ms. The reversed flow region near 
the lower wall has returned. The pressure gradient is again adverse. Re-examining figure 5.107, the 
regions of flow reversal roughly correspond to the regions where the flame is decelerating. When 
the flow reversal ahead of the flame is not present the flame is accelerating. The pressure gradient 
reversing is directly related to the flame motion. When the flame motion is decelerating, expansion 
waves are generated, causing an ad verse pressure gradient to form ahead of the flame. 1 

The centerline velocity at different times is shown in figure 5.116. The velocity at various channel 
locations ahead of the flame in the unburned gas is given in figure 5.117. Both figures show that the 
flow at the centerline of the channel ahead of the flame always remains positive, that is the flow is 
always directed toward the outflow of the channel. The centerline flow behind the flame is directed 
backwards toward the end wall, and can be greater in magnitude then the outward directed flow. 

Figure 5.118 gives the variation in the flame area in time. Shown for comparison is the flame 
area for case ADI. The flame area is defined for this case as the length of the 1000 K isotherm 
from the channel centerline to the closest point that it reaches from the lower wall. Comparing the 
variation of the flame area with the flame speed from figure 5.107, the area ratio follows the flame 


68 


speed. When the speed decreases, the area decrease, and vice-versa. 

The total chemical energy release and the total heat extracted by the cold wall are shown in 
figure 5.119. The total chemical energy release is defined as the integration of the local chemical 
energy release, as defined in equation (5.14), over the entire flame volume. Assuming a unit depth, 
the total chemical energy release can be calculated as 

E c = f J(-pq^)dxdy (5.16) 


This gives the total amount of energy at each time that is released in the chemical reaction. The 
total heat extracted from the flow is due to the cooling of the flow by the lower wall. Again assuming 
a unit depth, the total heat extracted from the flow is given by 


Qw = 



dy 


(5.17) 


where equation (2.33) has been used to evaluate the local heat flux to the lower wall. Due to the 
noncatalytic wall boundary condition, the term ( dY/dy) w is identically zero. Examining figure 5.119, 
the total chemical energy release follows the general shape of the flame area and speed variation. 
As flame area increases, there is more burning so the total energy release increases. The total heat 
removed by the cold wall is always increasing in time. This is due to the propagation of the flame 
constantly providing a longer heated zone behind the flame. 
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Figure 5.4: Temperature, density, pressure, and vorticity contours, 0.433 ms, ADI. 
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Figure 5.6: Temperature, density, pressure, and vorticity contours, 0.758 ms, ADI. 
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Figure 5.8: Temperature, density, pressure, and vorticity contours, 1.104 ms, ADI. 
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Figure 5.10: Temperature, density, pressure, and vorticity contours, 1.450 ms, ADI. 
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Figure 5.12: Boundary layer thickness for varying time, ADI. 



Figure 5.13: Displacement thickness for varying times, ADI. 
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Figure 5.14: Centerline u velocity profiles for varying times, ADI. 
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Figure 5.15: Boundary layer height for different x locations as a function of local time, ADI. 
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Figure 5.21: Temperature, density, pressure, and vorticity contours, 1.100 ms, AD2. 
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Figure 5.22: Velocity contours, velocity vectors and instantaneous streamlines, 1.100 ms, AD2. 
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Figure 5.23: Boundary layer thickness for varying time, AD2. 



Figure 5.24: Displacement thickness for varying times, AD2. 
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Figure 5.30: Flame position at tl 
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Figure 5.31: Flame velocity relative to laboratoi 
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Figure 5.32: Vertical position of the bulge in time, AD3. 



Figure 5.33: Length comparisons between flame points. 
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Figure 5.34: Temperature, density,pressure, and vorticity contours, 0.541 ms, AD3. 
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Figure 5.38: Temperature, density, pressure, and vorticity contours, 1.601 ms, AD3. 
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Figure 5.40: Temperature, density, pressure, and vorticity contours, 2.121 ms, AD3. 













Figure 5.45: Boundary layer height for different x locations as a function of local time, AD3. 



Figure 5.46: Boundary layer height as a function of the square root of the product of the local time 
and edge kinematic viscosity, AD3. 
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Figure 5.47: Flame position at the centerline and wall, AD4. 



Figure 5.48: Flame velocity relative to laboratory coordinates for the centerline and wall, AD4. 
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Figure 5.49: Temperature, density, pressure, and vorticity contours, 1.172 ms, AD4. 
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Figure 5.51: Displacement thickness at varying times, AD4. 



Figure 5.52: Boundary layer thickness at varying times, AD4. 
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Figure 5,57: Flame position at the centerline and wall, ADS. 



Figure 5.58: Flame velocity relative to laboratory coordinates at the centerline and wall. AD5. 
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Figure 5.59: Temperature, density, pressure, and vorticity contours, 1.284 ms, AD5. 
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Figure 5.63: Boundary layer thickness in local time, ADS. 
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Figure 5.64: Boundary layer thickness versus the square root of local time, ADS. 
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Figure 5.70: Position of the flame at the centerline and wall, ASl. 



Figure 5.71: Flame velocity relative to laboratory coordinates for the centerline and wall, ASl. 
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Figure 5.72: Relative flame area versus time, AS1. 
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Figure 5.73: Temperature, density, pressure, and vorticity contours, 0.102 ms, AS1. 
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Figure 5.74: Velocity contours, velocity vectors and 
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Figure 5.75: Temperature, density, pressure, and vorticity contours ,0.209 ms, AS1. 
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Figure 5.84: Boundary layer thickness for local times, ASl. 
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Figure 5.86: Position of the Same at the centerline and wall, AS2. 



Figure 5.87: Flame velocity relative to laboratory coordinates for the centerline and wall, AS2. 
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Figure 5.89: Temperature, density, pressure, and vorticity contours, 0.348 ms, AS2. 
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Figure 5.93: Temperature, density, pressure, and vorticity contours, 0.650 ms, AS2. 
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Figure 5.95: Temperature, density, pressure, and vorticity contours, 0.801 ms, AS2. 
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Figure 5.99: Boundary layer thickness at 
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Figure 5.103: Boundary layer growth comparisons for the spark ignition. 
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Figure 5.106: Flame position in time for the Isothermal wall case. 



Figure 5.107: Flame speed for the Isothermal wall case. 
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Figure 5.108: Temperature, fuel mass fraction, chemical energy release, and pressure contours 
0.217 ms, ISO. 
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Figure 5.110. Temperature, fuel mass fraction, chemical energy release, and pressure contours 
1.304 ms, ISO. 
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Figure 5.112: Temperature, fuel mass fraction, chemical energy release, and pressure contours 
1.760 ms, ISO. 
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Figure 5.114. Temperature, fuel mass fraction, chemical energy release, and pressure contours 
2.519 ms, ISO. 
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Figure 5.118: Time history of the relative flame area, ISO. 
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Figure 5.119: Chemical heat addition and thermal extraction in time, ISO. 
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Chapter 6 

Conclusions 


The interaction of a laminar flame with its self-generated boundary layer in a rectangular channel 
was numerically simulated using the two-dimensional, reacting, Navier-Stokes equations. A two 
species chemistry model was implemented which simulates the stoichiometric reaction of acetylene 
and air. Calculations were performed to investigate the effects of altering the boundary condition 
of the wall temperature, the Lewis number, the dynamic viscosity, and the ignition method. The 
purpose of this study was to examine the fundamental physics of the formation of the boundary 
layer and the interaction of the flame as it propagates into the boundary layer that its own motion 
has created. 

6.1 Summary of Results 

The pressure waves, that are emitted by the flame as it moves down the channel, induce a 
flow ahead of the flame in the same direction as the flame propagation. T his induced flow creates 
a boundary layer on the channel wall ahead of the flame. The flame then propagates into t his 
boundary layer that was created due to its own motion. These pressure waves are represented in the 
solution procedure. The compression and expansion waves produced by the flame are the mechanism 
by which the gas is set into motion. The interaction of the pressure waves with the gas is seen in the 
boundary layer growth ahead of the flame. Ahead of the initial pressure wave, which moves at the 
local speed of sound, the gas is unaware of the flame. The boundary layer begins to grow in tim e 
starting from the location of the initial pressure wave. This solution represents a direct numerical 
simulation of the acoustic interactions of the pressure waves created by the flame motion. 

The growth of the boundary layer is found to be self-similar in local time. Local time is defined 
to start at the instant the initial pressure wave generated by the flame reaches that point in the 
channel. The growth follows the general solution of Stokes’ First Problem. With the growth being 
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self-similar in local time, an analytical theory could be developed to calculate the flame position and 
speed in a channel. The boundary layer growth caused by the flame could be modeled as Stokes’ 
First Problem with an acceleration term. 

The deviation from the classical solution of Stokes’ First Problem is due to the flame acceleration. 
Flames that accelerate at a faster rate show more deviation then flames with a slower acceleration. 
Increasing the channel height has the effect of slowing the flame, thereby causing better agreement. 
When the viscosity was increased, the flame accelerated at a faster pace, causing a greater deviation. 
The reason for the acceleration of the flame is the upwelling at the wall. When the material burns in 
the flame, the density decreases. The material that burns near the wall is restricted by the boundary 
layer, which causes a flow blockage. The material that cannot expand ahead of the flame is forced 
to expand upwards, creating the upwelling. The upwelling turns at the centerline and provides 
additional momentum to the flame near the centerline. This creates the mushroom shape flame. 
When the channel height was increased, the upwelling was not dispersed over the entire flame, but 
acted on a flame region between the wall and the centerline. This created a rearward cusp at the 
centerline. Since the curved region near the wall grows in time, the upwelling increases, and its 
effects are dispersed over a greater portion of the flame. This leads to the rearward cusp at the 
centerline forming and then shrinking in time. 

When a simulated spark was used for the flame ignition, the flame propagated at a faster speed. 
This was due to the flame developing a greater flame front area in a shorter time as compared 
to the planar ignition cases. The greater the flame front area, relative to the channel height, the 
more material will be burned, and the more the flame will be accelerated. An interesting feature 
observed in the case with increased height was the formation of a flame on the end wall of the 
channel with a similar shape as the flame that propagates down the side wall. When the flame on 
the end wall reached the side wall a flame remnant formed that became detached from the main 
flame and eventually burned itself out. 

The results for an isothermal wall differed greatly from the results for an adiabatic wall. Since 
the wall temperature was fixed at the unburned gas temperature, the flame was quenched at the 
lower wall. Also, since the burned volume behind the flame was cooled by heat transfer to the 
lower wall, the burned material could expand in both directions. This had the effect of slowing the 
flame and at times even reversing the propagation direction. When the propagation direction was 
reversed, an area of reverse flow existed near the lower wadi. 
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6.2 Recommendations for Future Work 


The simulations presented in this study represent work done for a two-dimensional laminar flame 
in a rectangular channel in the absence of a gravitational force. Recommendations for future work 
include : 

• The gravitational term should be included. This would require that the orientation of the 
channel be specified, either vertical or horizontal. The flame shape and propagation will 
be altered by the buoyancy effects of the lower density burned gas behind the flame. In 
a horizontal channel, the buoyancy would require the simulation of the entire channel, the 
centerline symmetry condition could not be applied. 

• A calculation should be performed with a full-chemistry model, such as a methane-air or 
hydrogen-air model. A full chemistry model would allow the transport and chemical properties, 
such as the specific heats, to vary over the flame and ch ann el. 

• A turbulence model should be included. Turbulent flames propagate at a higher velocity 
then laminar flames. Also a turbulent boundary layer is thicker and would provide a greater 
blockage to the gas motion ahead of the flame. An alternate approach to adding an explicit 
model would be the direct numerical calculation of the turbulence by decreasing the grid size. 
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Appendix A 

Stokes’ First Problem with Various 
Accelerations 


Stokes’ First Problem was discussed in section 3.2.3 dealing with code validation. The general 
equation is derived from the x-momentum equation (2.2) with the following assumptions : 

• The u velocity is independent of the x direction, u = u(y,£). 

• v = 0. 


• The edge velocity is at most a function of time, U e = U e (t). 


With these assumptions, the x-momentum equation reduces to 

du 1 dp d 2 u ( A . 

■57 = -TT + (A* 1 ) 

dt pdx dy 2 

The edge velocity, U e (t), must satisfy the above equation. Substituting the edge velocity into 
equation A.l gives 

< ELl = (A. 2) 

dt p dx 

This equation can be used to replace the pressure gradient term with the acceleration of the edge 

flow, dU e /dt. The final equation to be solved is 

d 2 u 1 du ^ 1 dU c 

dy 2 v dt v dt 

with the boundary conditions given by 


(A.3) 


y = 0 : u = 0 

y — > 00 : u -4 U e (t) 


(A.4) 

(A.5) 


To show that flow with acceleration reduces the boundary layer thickness, three different accelera- 
tions will be discussed. 
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A.l Zero Acceleration 


This case has zero acceleration, giving a constant edge velocity, U e (t) = 
Stokes Fust Problem, as discussed in section 3 . 2 . 3 . Let 17 be defined as 

- U e - This problem is 

9- i 
2 yfM 

(A. 6 ) 

The solution to this problem is given by Schlichting as 55 as 


, , , x 

U' = 1 erfc(rj) 

(A.7) 

where erfc is the complementary error function defined as 


erfc(Tj) = 1 - erf ( 17 ) 

(A. 8 ) 

erf(»?) = - 7 = / exp(-r) 2 )dT] 
V* Jr, 

(A.9) 

This equation, along with the solutions for the next two sections, is shown in figure A.l. 

To obtain the boundary layer thickness, £ 99 , equation (A.7) is solved for the 77 where the velocity 
is 0.99% of the edge value. This gives the boundary layer thickness as 

699 = 3.64\/(l ft) 

(A. 10 ) 


A. 2 Constant Acceleration 


For this case the acceleration is a constant, giving a linear edge velocity profile. Let the acceler- 
ation be given by a, so that the edge velocity is given by U e (t) = at. The governing equation (A. 3 ), 


boundary conditions, and initial conditions are given as 

d 2 u 1 du _ a 

dy 2 v dt v (A ll) 

“(lb 0) = 0 (A. 12) 

U (M) = 0 (A. 13) 

u(y -¥ oo, £) -* at (A.14) 

To solve the problem, a change of variables is performed. 

v(y,t) = u(y,t) -at (A. 15 ) 
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(A. 16) 


Substituting into the above equations gives a new, homogeneous equation to be solved. 

dPv 1 dv _ 

dy 2 v dt 

with initial and boundary conditions given by 

v(y,Q) = 0 (A. 17) 

u(0,0 = - at (A.18) 

v(y -4 oo, <) — ► 0 (A. 19) 

A general solution exists for the problem with the boundary condition at the wall given by 

v(0,t) = K<f>^ n ,n — —1,0, 1,2,... (A.20) 

The solution is given by Carslaw and Yeager 69 as 

v = AT(in + l)(4t)* n i n erfcfa) (A.21) 

where T is the gamma function. To complete the solution, the following formula are used to obtain 
a solution without the complex variable, i. 


2m n erfc(x) = 

t n 2 erfc(x) - 2xx n ! erfc(x) 

(A.22) 

x 2 erfc(x) = 

^ [erfc(x) - 2xxerfc(x)] 

(A.23) 

ierfc(x) = 

—= exp(~x 2 ) — xerfc(x) 
V* 

(A.24) 


For the problem as stated in equations (A. 16 - A. 17), K = —a and n = 2. The gamma function 
for n = 2 is given by T(2) = 1. Solving the transformed problem and transforming back to the 
original variable, u(y,t), gives 

— = 1 - (1 + 2r? 2 )erfc (rj) + -^=rj exp(-T) 2 ) (A.25) 

at y/7T 

This equation is shown in figure A.l. To obtain the boundary layer thickness, this equation is solved 
for the 7} where u(y, t)/at = 0.99. This value is r/ = 1.45. The boundary layer thickness for the 
constant acceleration problem is given by 

£99 = 2.90v^i t = (0. 797)3.64 (A.26) 

T his shows that for a constant acceleration, the boundary layer thickness grows slower then for the 
constant velocity case. 
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A. 3 Linear Acceleration 


For this case, the velocity follows a parabolic distribution which gives a linear acceleration. Let 
the acceleration be given by at, so that the edge velocity is given by U e (t) = 0.5 at 2 . The governing 
equation (A.3), boundary conditions, and initial conditions are given as 

dPu \du at 

dy 2 ~udi~~7 (A ' 27 ) 


u(y,0) = 0 

u(0,t) = 0 

u(lf-Kx>,t) hat 2 

As was done in the previous section, a change of variables is performed. 

v(y,t) = u( j/.O-io * 2 

Substituting into the above equations gives a new, homogeneous equation to be solved. 

d 2 v _ 1 dv _ 
dy 2 v dt 

with initial and boundary conditions given by 

v(y, 0) = 0 
v(0,t) = -hat 2 

v(y oc,t) 0 


(A.28) 

(A.29) 

(A.30) 


(A.31) 


(A.32) 


(A.33) 

(A.34) 

(A.35) 


Making use of the general solution given in equation (A.21), K = -\a, n = A, and the gamma 
function for n — 4 is given by T(3) = 2. Solving the transformed problem and transforming back to 
the original variable, u(y, t), gives 


T^T = 1 ~ + 2»7 2 )erfc(ij) - -h=r)exp(-T} 2 ) 

+ f* ex P(~V) ~ »?erfc(»;)J 
~ 3^ + 2r? 2 )erfc(jj) - -^Ljjexpf-j? 2 ) 


(A.36) 


This equation is shown in figure A.l. To obtain the boundary layer thickness, this equation is solved 
for the t) where u(y,t)/U e = 0.99. This value is t) = 1.23. The boundary layer thickness for the 


constant acceleration problem is given by 


<$99 = 2.46v^ = (0.677)3.64v^t 


(A.37) 
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This shows that for the case of a linearly accelerating edge velocity, the boundary layer grows slower 
then both the constant and zero acceleration cases. 



Figure A.l: Velocity profiles with varying edge accelerations 
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Appendix B 

Density and Temperature 
Variations Due to a Compression 
Wave 


A finite compression wave is defined as a wave where the changes, or perturbations, from the 
initial conditions may not be small. This is opposed to a weak or acoustic compression wave where 
the perturbations are assumed to be small so that linear relationships may be developed. When a 
finite compression wave passes through a point in space, the pressure, density, and temperature all 
increase. Let the initial condition be given by pi, Tj, and pi. After the wave has passed by, the new 
conditions in region 2 are given by a perturbation from the initial conditions. 


P2 = 
T 2 = 

P7 = 

Assume the wave is isentropic. This allows t] 
isentropic flow relations. 


Pi + Ap 

(B.l) 

7i + AT 

(B.2) 

Pi + A p 

(B.3) 


conditions in regions 1 and 2 to be related by the 


tey-GT* 


(B.4) 


where is the ratio of specific heats. Substituting equations (B.2) and (B.3) and redu cin g gives 

/ Ap \ ( AT 

l^ +1 ) = hr +1 ) (B - 5 > 

To relate this result to the relative pressure increase, the equation of state is used. In regions 1 
and 2, the equation of state becomes 


Pi = piR.Ti 


(B.6) 
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P2 = P2-R«?2 

Using equations (B.l - B.3), the equation in state 2 becomes 

Pi + Ap = (pi 4- Ap)R,(Ti + AT) 

Pi + A p — piR,Ti + pi R/AT + ApR t Ti + ApR t AT 

Dividing both sides by equation (B.6) and reducing using equation (B.6) gives 

Ap _ Ap AT A p AT 

Pi Pi Ti pi Ti 


(B.7) 


(B.8) 

(B-9) 


(B.10) 


This equation provides the relationship between the relative increases of the density and temperature 
as related to the pressure increase. 

B.l Acoustic Wave 


Assume that the ratio of the perturbation to the initial condition is small. 

AT 


Ti 


« 1 


(B.ll) 


(B.12) 


Using the binomial expansion, neglecting higher order terms, equation (B.5) reduces to 
For the present chemistry model, 7 = 1.25 so that the above equation becomes 

7T i4 (f) <B13) 

This shows that for a finite compression wave, the relative increase in density is greater then the 
relative increase in temperature, 4x for the present chemistry model. Assuming that the product 
of perturbations is small, equation (B.10) becomes 


Ap . A p AT 
Pi ~ Pi + Ti 

(B.14) 

Ap . 7 AT 

Pi ~ 7 - 1 ?i 

(B.15) 

Ap . AT 
Pi 5 Ti 

(B.16) 
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B.2 Finite Wave 


For a finite wave the perturbations need not be small. Solving equation (B.5 for the density 
variation gives 


a p ( at \ 1/(7-1) 

Pi~\T x +1 ) -1 

Substituting this relation into the pressure equation gives 

A p (AT \ ■*/('»-!> 

Pi “ (ri +1 ) 


(B.17) 


(B.18) 


For a given relative increase in the pressure due to a finite wave, the relative increase in the density 
is greater the the relative increase in the temperature. 
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Appendix C 

Laminar Flame Theory 


C.l Overview 

Laminar flame theory provides a means to estimate the flame speed, which is defined as the speed 
of the flame relative to the flow ahead of it, and the flame thickness. Mallard and LeChatelier 27 
developed a model of a flame as shown in Figure C.l. The flame is divided into two region. The first 
region, Zone I in the figure, is a preheat region. There is no chemical reaction in this region. The 
unburned material is heated due to thermal conduction up to a temperature T*, which is assumed to 
be close to the temperature of the burned material, 7V In the second region, Zone II in the figure, 
the material combusts. The flame is stationary with the gas flow coming towards it at the laminar 
flame speed, 5/. 

C.2 Laminar Flame Speed 

Zeldovich and Frank-Kamenetskii 60 " 62 and Semenov 63 developed a laminar flame theory based 
on the model described above. They solved the species mass conservation equation and the energy 
conservation equation to determine a formula for the laminar flame speed. The assumptions that 
they made are : 

• One-dimensional (~ = 0). 

• Steady state (m = °)- 

• Viscous stress terms are negligible. 
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Figure C.l: Temperature profile across a laminar flame 

• Constant pressure. 

• The number of moles is constant during the reaction. 

• The specific heat, Cp, is constant 

• The thermal conductivity, is constant 

• Lewis number is unity. From equation (2.30) this implies that k/cp = pD. 

With the above assumptions, the global continuity equation, (2.1), can be written as 

dpu 

— 0 => pu = constant (C.l) 

This is a statement that the mass flow rate is constant. Applying this equation to the flame model 
in Figure C.l gives 


pu = p u Si 


(C.2) 


For a steady, one-dimensional flame, the species continuity equation, (2.5), written in terms of 
the fuel mass fraction, Y , is given by 


pu 


dY 

dx 


dx 


, dY x 


4- u) 


(C.3) 
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where u is the production rate of fuel. Using equations (2.18) and (2.19), this terms becomes 

u = p 2 YAexp(-T a /T) (C.4) 


Transforming this equation using <p = Y — T u , equation (C.3) can be written as 

n cP<p dtp 


(C.5) 


The energy equation, (2.4), is written in terms of temperature. For the current chemistry model, 
this equation becomes 


,k x dPT dTuq 
Cp dx 2 ^ dx Op ~ 

Transforming by 6 = (cp/q)(T u - T), this equations becomes 

,k,(Pe dB 

(-)ti ~ + v = 0 

Cp dx 2 dx 

The boundary conditions for equations (C.5) and (C.7) are given by 

1 = -»{ 8 = 1 

{.*: 


X = +00 


' = 0 
0 

= ~Y U = -1 
(cp/q)(T u - T b ) 


(C.6) 

(C.7) 

(C.8) 

(C.9) 


Comparing equations (C.5) and (C.7) and the above boundary conditions, it can be seen that the 
solutions of the two equations will be the same, 4> = 0, if the boundary conditions at infinity are the 
same. 


-Yu = ^(T u -T b ) 
Q 

CpTb = CpT u -f q 


(C.10) 


This condition is satisfied if the flame is adiabatic. To prove this, examine equation (2.33). Setting 
the x component of heat flux to 0, and integrating over the flame from T u to T, we get 

„ t dT n dY 


-k f dT = pDq f Y dY 
J T u J K 

CpT + qY = CpT u + q 


(C.ll) 


The left hand side of this equation is evaluated at the burned side, where T = T b , and Y b = 0, to 
get 


OpT b — CpT u + q 


(C.12) 
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This equation is identical to equation (C.10). Making the assumption that the flame is adiabatic, 
equation (C.ll) can be used in place of equation (C.5). This reduces the system to one differential 
equation and one algebraic equation. 

The solution procedure is to integrate equation (C.6) over the two zones in figure C.l and match 
the first derivative at x = 0. This ensures that the heat flux is balanced between the two regions. 
In Zone I there is assumed no reaction, lj = 0. 


#T 

(PUCp) 


dx 2 

V k ) 

1 dx 


(C.13) 


with boundary conditions for this region 


r - 00 , T = r„, f =0 
10- T = 


Ti 


(C.14) 


Integration this equation from — oo to 0 and applying the boundary conditions at — oo gives 


— = (mCp (T — T ) 
dx k U lu) 


(C.15) 


This equation for the derivative is evaluated at x = 0~. It is assumed that T, is very close to T b , so 
that Ti « T;,. 


(CL-= ( ci6 > 

In Zone II, it is assumed that the convective term in equation (C.6) is negligible, since Ti is very 
close to 2V 


<PT uq 
dx 2 + k 


= 0 


(C.17) 


with boundary conditions for this region given by 


f oo, T = Ti, 
\ 0+, T = T, 


f =0 


(C.18) 


Integrating the above equation from 0 + to oo, applying the boundary conditions, and evaluating the 
result at 0 + gives 



(C.19) 


where 


1 = 



(C.20) 


The lower integration bound cam be replaced by T u since it is assumed that there is no reaction from 
Ti to T u . 
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At x = 0 the two derivative equations, (C.16) and (C.19), must be equal to ensure a heat flux 
balance between the two regions. 

(C.21) 

Solving the above equation for pu, substituting the result into equation (C.2), and reducing using 
equation (C.12), gives an equation for the laminar flame speed, Si as 

Si = 

According to Kuo, 2 this formula doesn’t give very accurate results, but it does predict the trends 
well. One of these trends is that the flame speed is proportional to the square root of the thermal 
conductivity, Si a y/k. 

C.3 Laminar Flame Thickness 




Turns 64 outlines a method from Spalding 65 to calculate the laminar flame thickness, defined by 


ft -ft 


(C.23) 


Spalding 65 assumed that the temperature gradient through the flame was constant, making the 
temperature profile a straight line. Equation (C.6), the energy equation, is integrated from -oo to 
+oo, with the following boundary conditions 


f +oo, T = ft, 

\ -oo, r = r u , 


dT/dx = 0 
dT/dx = 0 


This integration yields 


-«-i£ 


(C.24) 


(C.25) 


Solving equation (C.23) for dx , substituting into the above equation, reducing using equations (C.12) 
and (C.20), and substituting in for pu from equation (C.2), gives 


p u Si(Tt - T u ) = 6 t I 


(C.26) 


The term I can be removed from this equation by solving equation (C.22) for I and substituting. 
After reduction, the equation for the laminar flame thickness becomes 


= *(-*-) 
Si VPtiCp/ 


(C.27) 


Since the laminar flame speed was shown in the previous section to be proportional to the square 
root of the thermal conductivity, this equation shows that the thickness is also proportional to the 
square root of the thermal conductivity, St a y/k. 
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C.4 Quenching Distance 


Assume that the flame propagates between two parallel plates separated by some distance. 
The quenching distance, d q , is defined as the smallest plate separation that will allow the flame 
to propagate. If the plate distance is set any smaller, the flame will be quenched. The quenching 
distance can be calculated by equating the heat generated due to the chemical reactions, q c , to the 
heat lost at the wall, q w , 2 ’ 64 

Let A w represent an area of wall surface. The total heat loss to both walls in given by 

Qw = 2k w (C.28) 

Let the temperature profile is linear from the wall to the center of the channel. Also assume the the 
wall temperature, T w is equal to the unburned gas temperature, T u . Using the linear temperature 
profile, the heat loss at the wall can be written as 


Qw = 2k w A w i— “ (C.29) 

The heat release due to the chemistry is estimated as the average reaction rate over the volume 
between the plates, times the heat of formation. The reaction rate, 7, is given by equation (C.20) 

Qc = jT ^ T " Q A w d q (C.30) 


Equations (C.29) and (C.30) are equated and solved for the quenching distance, d q . 

j2 4^m(T; T u )(7 , — Tu) 

Iq 


(C.31) 


Solving equation (C.22) for 7, and making use of equation (C.27), the reaction rate is given by 


r _ 2 k(T b - T u ) 
~ erf 


(C.32) 


Substituting this result into equation (C.31), and using equation (C.10) for the heat of reaction, q, 
the quenching distance is calculated as 





(C.33) 


In this equation, T q is the lowest temperature at which the flame can propagate. Also k is evaluate 
at the same temperature used in the determination of the laminar flame speed. 
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Appendix D 

General Outline of CMRFAST2D 


This Appendix will outline the basics of CMRFAST2D, a Connection Machine, Reacting, FAST, 
Two-Dimensional code. A flowchart of the code structure is shown in figure D.l. 

(i) The initial values of the flowfield properties are defined either from a restart file or from user- 

defined initial values. The boundary conditions are established and the gridding is defined. 
The viscous stress terms, heat conduction terms and mass diffusion terms are initialized to 
zero. 

(ii) The timestep is determined as discussed in section 2.5 The total time of the simulation is 

calculated. 

(iii) The viscous stress terms in both the momentum and energy equations are calculated using a 
second order method. 

(iv) The heat conduction terms in the energy equation are calculated using a second order method. 

(v) The mass diffusion terms in the species continuity equations and the heat flux due to the mass 

diffusion in the energy equation are calculated using a second order method. 

(vi) The governing equations are integrated in time using the LCPFCT 40 routine as described in 
section 2.4. The viscous, heat conduction, and mass diffusion terms are used as source terms. 
The integration technique uses directional splitting to integrate in the x and the y coordinate 
directions. 

(vii) The chemical rate equations are integrated in time, as described in section 2.6. 

(viii) This routine updates the values at the boundaries using the characteristic boundary method 
described in section 2.7. Boundaries can be specified as subsonic outflow, inflow, solid wall, 
symmetry, or extrapolation (supersonic outflow). 
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(ix) This routine outputs the flowfield values at selected times. It also generates a restart file at 
selected times. 

(x) If the total time is less then a stop time, the loop continues by deter minin g a new time-step. If 

the time is greater, a restart file is written and the program ends. This loop criteria can also 
be based upon an iteration counter. 
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(i) 

(ii) 

(Hi) 

(iv) 

(v) 

(vi) 

(vii) 
( viii ) 
(xi) 
(x) 


Figure D.l: CMRFAST2D program flowchart. 
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